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COLD-SAT Dynamic Model 

Neil S. Adams 
Analex Corporation 
Brook Park, Ohio 44142 

and 


Gary Bollenbacher 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


1.0 Introduction 


The COLD-SAT dynamic model is a six-degree-of-freedom, rigid-body computer simulation of a 
spacecraft in orbit around the Earth. The model was developed as part of a conceptual design review 
at % the NASA Lewis Research Center to study the feasibility of a proposed spacecraft known as 
COLD-SAT. COLD-SAT is an acronym for cryogenic on-orbit liquid depot— storage, acquisition, 
and transfer. The purpose of the COLD-SAT spacecraft is the investigation of the fluid dynamics and 
thermodynamics of subcritical cryogenic fluids in a microgravity environment. In order to carry out 
these investigations, a series of experiments will be performed onboard the spacecraft during its 
proposed 60-day operational life. 


As envisioned during the study, the COLD-SAT spacecraft will be placed in a circular low Earth 
orbit at an initial altitude of 550 nautical miles and an inclination of 18°. The spacecraft will 
incorporate three hydrogen tanks: one large tank (the supply tank), which will be filled with liquid 
hydrogen at lift-off, and two smaller tanks (the large receiver tank and the small receiver tank), whic 
will be filled with and drained of liquid hydrogen during the course of various onboard experiments. 

In order to provide a controlled microgravity environment during the execution of some of the 
experiments, axial thrusters that provide three levels of microgravity acceleration will be mouted on 
the spacecraft. The cumulative time of this thrust over the 60-day mission will be approximately 
60 hours; for the remaining time the spacecraft will be in a zero-gravity condition. During the course 
of some onboard experiments unwanted disturbances (resulting primarily from attitude control systern 
activity) must be minimized. Otherwise these disturbances could be transferred to the cryogenic fluids 
and degrade the quality of the resulting experimental data. This effect is of greatest concern during 
experiments involving the application of axial thrust. The design life of the spacecraft was specified 

as six months. 


At the outset of the design study the baseline attitude control system was an all-thruster system 
consisting of six pairs of hydrazine thrusters. A single-axis reaction wheel was explored for meeting 
the requirements for low microgravity disturbances. A two-axis gimbal of one axial thruster was also 
provided to augment the fixed thrusters. As the study progressed, the paired control thrusters on two 
of the spacecraft axes were replaced with single thrusters to reduce system cost and complexity. 


The design of the COLD-SAT dynamic model (fig. 1.1) was shaped largely by the requirements of 
the COLD-SAT design study. The model was a major tool that was used throughout the study, 
primarily to evaluate and compare the performance of several attitude control system design options. 
Other areas of investigation for which the model was extensively employed were the following: 



(1) The analysis of the microgravity environment within the liquid hydrogen tanks, which was 
largely produced by the activity of the attitude control system 

(2) The evaluation of the effects of prolonged periods of constant, low-level thrust on the 
spacecraft orbit over the 60-day mission life 

(3) An assessment of the effects of liquid hydrogen sloshing in each of the three tanks on the 
spacecraft dynamic behavior 

Although throughout its development the design of the COLD-SAT dynamic model was guided by the 
requirements of the design study, an attempt was made to keep the model as general as possible so 
that its usefulness would extend to future projects as well as to COLD-SAT. To achieve this goal, 
approximations and simplifications were avoided to the maximum extent feasible. 

The resultant model consists of three parts: a translation model, a rotation model, and a slosh 
model (figs. 1.2 to 1.4). Each of the three parts is described in detail in subsequent sections of this 
report. The translation model simulates the motion of the spacecraft center of gravity about the Earth 
under the influence of gravitational force, atmospheric drag, and the thrust produced by the axial 
thrusters and by uncoupled control thrusters. The rotation model simulates the attitudinal motion of 
the spacecraft about its center of gravity under the influence of various disturbance torques and control 
torques acting on the spacecraft. The slosh model computes the torque on the spacecraft produced by 
the motion of the liquid hydrogen in any or all of the three tanks. In addition to the primary function 
of each part of the model, other computations are included as required. Thus, the translation model 
also computes the desired spacecraft attitude for a large number of attitude options, and it computes 
the torque resulting from axial thrust misalignment acting on the spacecraft. The rotation model 
simulates the attitude control system and computes the microgravity environment within the spacecraft 
body. 

Some of the features included in the model are 

(1) The capability to simulate a large number of possible spacecraft attitudes 

(2) An Earth gravity model that includes the fundamental component and zonal harmonics up to the 
fourth and accurately represents nodal regression and apsidal progression 

(3) Accurate simulation of gravity-gradient torque and atmospheric density 

(4) A tilted-dipole model of the Earth’s magnetic field, which provides sufficient accuracy over 
the range of orbital altitudes used for COLD-SAT 

(5) A Sun model that gives the position of the Sun as a function of time and date and determines 
if the Sun is blocked by the Earth 

(6) A two-axis slosh model that can simulate slosh in one tank or in more than one tank 
simultaneously 

(7) An axial thrust model that includes thrust decay due to the unregulated (blowdown) system 

(8) Simulation of control thnister (reaction control system, or RCS) startup and shutdown 
transients, as well as minimum on-time constraints 


2 


(9) A reaction wheel model that includes viscous and static friction 

(10) A microgravity acceleration computation that includes gravity-gtadient-induced acceleration, 
as well as accelerations caused by spacecraft lateral and rotational motion and angular 
acceleration 

(1 1) Quaternions that are used throughout the model to represent spacecraft attitude, which result in 
a representation that is valid for all rotations, is devoid of singularities and small-angle 
approximations, and requires a minimum number of numerical integrations per time step 

Other computations that are included in the model but are modeled in a simplified or approximated 
form or are only modeled for certain attitude options are 

(1) Solar pressure torque, which is approximated so that it is only valid for the attitude selected 
for the COLD-SAT spacecraft with small attitude errors 

(2) Aerodynamic drag torque, the accuracy of which is limited by the modeling of the 
spacecraft surfaces as simple geometric shapes 

(3) Magnetic torque, the accuracy of which is limited by the modeling of the Earth’s magnetic 
field as a simple dipole 

However, these torques are very small relative to the dominant disturbance torques so that accurate 
simulation was not required. Furthermore, highly accurate computation of these torques would result 
in a model that required excessive computer run time. 

Another limitation of the model is the computation of slosh torque, which cannot be simulatedfor 
certain attitude options because the slosh model in its present form will allow no more than a 180^ 
spacecraft rotation with respect to inertial space. Thus, the slosh model will only function properly 
with attitude options that are nearly inertially fixed. 

Often during the design phase of the COLD-SAT study, alternative systems were evaluated to 
determine which system best met COLD-SAT’s requirements in the most cost-effective manner. 
Consequently, the COLD-SAT dynamic model was constructed to have the capability to simulate all 
possible alternative systems. Examples of this include the following items. 

(1) Both regulated and blowdown thrust systems were considered for COLD-SAT, and both can be 
simulated by using the COLD-SAT dynamic model. 

(2) A single-axis reaction wheel was considered during the design phase to improve the 
microgravity environment during periods of no axial thrust. Although the reaction wheel was 
eliminated from the final design, it is included in the model. 

(3) Both coupled and uncoupled RCS thruster systems were considered for COLD-SAT. An 
uncoupled scheme was chosen for COLD-SAT to satisfy cost and complexity considerations, but the 
model can simulate either system. 

The model is implemented by using the EASY5 modeling and analysis package and will run on a 
variety of computing platforms, including the Cray X/MP supercomputer. Execution time is a strong 
function of the thrust level simulated. Simulation of the full-up model for one complete orbit on the 
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Cray, using the most efficient numerical integration method available, varies from about 1 min for no 
thrust to about 12 to 15 min for the highest thrust level. Execution time for the translation model 
alone, simulating the entire 60-day mission, is on the order of 2 to 3 hours of Cray time. 

The remainder of this report provides detailed descriptions of all aspects of the model. Included are 
the mathematical formulations and definitions of variables. The constants, quaternion algebra, and 
coordinate systems used are discussed in the appropriate appendixes. The code listings are available 
through the Computer Software Management and Information Center, Athens, Georgia, 30602. 



Figure 1.1 .—Overview of COLD-SAT model. 
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Figure 1 .4.— Slosh model. 
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2.0 Translation Model 


The translation model propagates the position and velocity vectors from an assumed initial on-orbit 
state for a spacecraft under the influence of various forces. The initial state is computed in the model 
as a function of the following user-specified orbital elements: 

(1) Perigee altitude 

(2) Eccentricity 

(3) Inclination 

(4) Argument of perigee 

(5) Right ascension of the ascending node 

(6) Time since perigee passage 

Other inputs to the translation model are the initial date and time (in GMT) used to compute the 
position of the Sun and the desired attitude option and orientation. A block diagram of the translation 
model is shown in figure 1.2. 

Three rectangular, right-handed coordinate systems are used for the translation model. They are 
discussed briefly here and in more detail in appendix C: 

(1) True-of-date (TOD) coordinate system (fig. 2.1): A pseudo-inertial, equatorial, Earth-centered 
system whose X axis points in the direction of the first point of Aries (the vernal equinox) and whose 
Z axis points in the direction of the North Pole. The axes of this system are designated X, Y, and Z 
See figure C.l. 

(2) Spacecraft-body system: A coordinate system fixed with the center-of-gravity point (CG) of the 
spacecraft, where the origin is the actual CG location of the spacecraft. The x axis is defined as 
parallel to the long axis of the spacecraft, pointing toward the fore end, the y axis is parallel to the 
solar panel connecting struts, and the z axis points in the direction opposite the spacecraft high-gain 
antenna. The axes of this system are designated x, y, and z. See figure C.2. 

(3) Alternative spacecraft-body system: A coordinate system similar to (2) with the origin located 
where the plane containing the rear thrusters intersects the line that is parallel to the long axis of the 
spacecraft and passes through the nominal CG point (thus, if the y and z CG offsets from the nominal 
are zero, then y = y' and z = z’). The axes of this system are specified as x\ y' , and zf. See 

figure C.2. 

All forces acting on the spacecraft, excluding gravity, are computed in spacecraft-body coordinates and 
then transformed into TOD coordinates. The gravitational force is computed directly in TOD 
coordinates, and the integration of the resulting accelerations from all applied forces is carried out in 
the TOD coordinate system. 

Because the model was required to accurately simulate the spacecraft orbit for periods of time 
extending up to 60 days, changes in the spacecraft mass were also modeled. The mass changes of the 
COLD-SAT spacecraft are significantly larger than those associated with a typical spacecraft because a 
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high percentage of the spacecraft mass is composed of expendable fluids, which include liquid 
hydrogen in addition to hydrazine (N 2 H 4 ). 

The translation model includes a number of other auxiliary functions that are not directly associated 
with the computation of the spacecraft position and the velocity vectors. These include the 
computation of the desired attitude matrix and subsequent extraction of an equivalent quaternion and 
the computation of the Sun vector and a vector normal to the ecliptic plane. These two vectors are 
required for the computation of a desired attitude quaternion for certain attitude options. Also 
included in the translation model are computations of the Greenwich hour angle (GHA) and the 
conversion of the calendar date and the Greenwich mean time (GMT) to the Julian date. 

The remainder of this section provides detailed descriptions of the various 
elements of the translation model as follows: 

(1) Conversion of initial orbital elements into an initial state vector in TOD coordinates 

(2) Computation of the Julian date, the Sun vector, the ecliptic-normal vector, and the GHA 

(3) Computation of all simulated forces acting on the spacecraft and transformation of these 
forces to TOD coordinates by using the attitude quaternion 

(4) Approximation of the spacecraft mass 

(5) Integration of the resulting accelerations to obtain velocity and position vectors 

( 6 ) Computation of the desired spacecraft attitude quaternion 

(7) Computation of torque due to axial thruster misalignment 

( 8 ) Computation of classical or osculating orbital elements from the state vector for printing and 
plotting 


2.1 Interfaces of Translation Model 

Interfaces between the translation model and other parts of the overall COLD-SAT dynamic model, 
as well as input and output quantities, are as follows: 

( 1 ) The following quantities are computed by the translation model and passed to the rotation 
model: 

(a) The position vector is used for disturbance-torque computations. 

(b) The atmospheric drag forces are used for (aerodynamic drag) disturbance-torque 
computations. 

(c) The Sun position is used for (solar pressure) disturbance-torque computations. 

(d) The desired attitude quaternion is used to compute attitude error. 
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(e) The GHA is used for (magnetic) disturbance-torque computations. 

(f) The desired attitude quaternion at the beginning of the simulation is used to initialize the 
actual spacecraft attitude quaternion. 

(2) The desired attitude quaternion at the beginning of the simulation is passed to the slosh model 
to provide a rotation reference. There are no other interfaces between the translation model and the 
slosh model. 

(3) The following quantities are computed by the rotation model and passed to the translation 
model: 

(a) The actual spacecraft attitude quaternion, which is used to transform forces acting on the 
spacecraft from body coordinates to TOD coordinates. (If the translation model is used alone, the 
desired attitude quaternion is used in place of the actual attitude quaternion). 

(b) The translational RCS accelerations resulting from unpaired thrusters. (A constant, average 
value is used to represent the translational RCS force resulting from uncoupled thrusters if the 
translation model is used alone.) 

(4) The following quantities, which are specified by the user at the start of simulation, are input to 
the translation model: 

(a) Initial orbital elements (see beginning of section 2.0) 

(b) Initial calendar date and GMT 

(c) Spacecraft attitude option and axis alignment (see section 2.6) 

(5) The primary outputs of the translation model are numerous classical or osculating orbital 
elements, which can be plotted or presented in tabular form as a function of time. 


2.2 Computation of Initial Position and Velocity Vectors From Given 
Orbital Elements 

An initial state vector, which consists of the initial position vector Rq and the initial velocity vector 
v 0 , must be specified at the beginning of the simulation. It is computed from the specified orbital 
elements. This section is concerned with this computation. The following standard orbital elements 
are used as initial conditions. From these values the initial position and velocity vectors can be 
computed. Some of these orbital elements are illustrated in figure 2.2. 

a semimajor axis, ft 

e eccentricity 

i orbit inclination from equatorial plane, deg 

Q right ascension of ascending node, deg 
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co argument of perigee, deg 


T p time since perigee passage, sec 


From these given elements, additional orbital elements are computed as follows: The mean anomaly 
(in radians) is given by 


M = 



where Ji is the Earth gravitational constant (see appendix A). The eccentric anomaly E may be 
computed by the following iterative pseudocode: 


E = M + e sin M 
Loop: M e = E - e sin E 

If (abs(M - M e ) < TOL) STOP, answer is E. 
E = E + (M - M e )/(l - e cos E) 

Go to loop 


where TOL is a small number (TOL « 1) that indicates the accuracy of the result. 


An Earth-centered, right-handed coordinate system is defined with x and y axes lying in the 
spacecraft orbit plane. The x axis is in the direction of orbit perigee, and the z axis is in the direction 
of the spacecraft angular momentum vector. This coordinate system is called the spacecraft orbit 
(SCO) system. Refer to appendix C for a detailed description of all coordinate systematised in the 
model. The position and velocity vectors in spacecraft orbital coordinates, R sco and v sco , 
respectively, are computed as follows: 


R, 


SCO 


a(cos E - e) 


a sin £yl - e 2 
0 


sj\ia sin E 
R 

\[^a( 1 - e 2 ) cos E 
R 

0 

Position and velocity vectors are then transformed to TOD coordinates as follows: 



Ro = [SCO - TOD]/? sco 
7 0 = [SCO - TOD]7 sco 
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where 


[SCO - TOD] = 


cos O sin Q 0 
-sin O cos Cl 0 
0 0 1 


1 0 0 
0 cos i sin ( 
0 -sin i cos i 


cos (o sin to 0 
-sin co cos to 0 
0 0 1 


and 


Q right ascension of ascending node, deg 
i inclination, deg 
co argument of perigee, deg 


There are two special cases of spacecraft orbit where one or more of the input orbital elements are 
not defined: 


Case 1: e = 0, / * 0: If the spacecraft orbit is circular, orbit perigee is undefined, and the 
orientation of the SCO x axis is arbitrary. For a circular orbit the argument of perigee input co 
should be set to zero, so that the SCO x axis is coincident with the ascending node. Thus, the T p 
input is interpreted, in this case, as the time since the spacecraft passed the ascending node. 

Case 2: e = 0, i = 0: If the orbit is circular and lies in the equatorial plane, both the argument of 
perigee and the right ascension of the ascending node are undefined. In this case, both of the input 
quantities (co,0) should be set to zero, so that the SCO coordinate system becomes identical to the 
TOD system and the T p input is interpreted, in this case, as the time since the spacecraft passed the 
TOD X axis. 


2.3 Computation of Forces Acting on Spacecraft 


13.1 Gravity 

The Earth gravity model computes the linear acceleration of the spacecraft due to the gravitational 
field of the Earth. The fundamental component (spherical-Earth component) and as many as three 
additional terms resulting from the oblateness of the Earth (zonal harmonics) can be computed, so that 
nodal regression and apsidal progression can be accurately simulated. 

The gravitational potential function of the Earth (from ref. 1) is as follows: 

2 3 4 

J + ^Earth (1 _ 3 9in 2 e) + ^Earth (3 _ 5 sin 2 ©) sin 0 - (3-30 sin^ + SS 
2R 2 2 R 3 8 R 4 

where 


0 >(/?, 0 ) = 


R magnitude of spacecraft position vector, ft 
6 geocentric latitude of spacecraft 
Earth Earth gravitational constant, 1.4076469 * 10 16 frVsec 2 



^Eaiih Ea^ equatorial radius, 20 925 656 ft 
J 2 second zonal-harmonic coefficient, 1.0827 * 10 
J 3 third zonal-harmonic coefficient, -2.56 x 10 6 
fourth zonal-harmonic coefficient, -1.58 x 10 6 

The terms J v and J 4 in the preceding equation are called zonal harmonics and result from the 
oblateness (flattening of the poles) of the Earth, which is the cause of the nodal regression of the orbit 
plane and the apsidal progression of the argument of perigee. Tesseral harmonics (resulting from 
longitudinal variations in the density of the Earth), zonal harmonics beyond the fourth, and solar and 
lunar gravitational effects were considered negligible for these purposes and are not modeled. 


The gravitational acceleration can be found by taking the gradient of the potential function as 
follows: 


V<t>(/?,0) 


- d<t>(R, 6 ) , -» 1 

,0) = \> R g R + uege = v R — - + «e - 


1 0O(/?,0) 


00 


where \? R and tT e are unit vectors in the R and 0 directions in spherical coordinates, respectively. 
Refer to figure 2.1 for an illustration of this coordinate system. 

By evaluating the appropriate partial derivatives, the radial and latitudinal components of the 
gravitational acceleration (in spherical coordinates) can be determined as follows: 

- - i!52H - 2 - 3 sin2e ) - ~ 5 0 - 5 sin 0 

♦ 4 - 30 sin2e * 35 5i " 40) 

O 


g e = 0 cos e ^ CO, e - 15 sin 2 !) cos 6) 

- 4 e COS 6 * 140 ,in 3 0 co, 6) 

o 

As there is no longitudinal variation in the potential function, the longitudinal component of 
gravitational acceleration is zero. 

S<j> = 0 

The gravitational acceleration of the spacecraft in spherical coordinates is then converted to 
Cartesian TOD coordinates as follows: 


Q 


cos r| cos 0 -cos q sin 0 

r - 

OX 



8r 

8y 

= 

sin q cos 0 -sin q sin 0 

£e 



sin 0 cos 0 



where q is the right ascension of the spacecraft position vector. 
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Right ascension and latitude are computed from the spacecraft position vector /? as follows: 


0 = tan 


n - 

(Note that all arctangent functions used in this document are of the four-quadrant form, which 
considers the numerator and denominator of the argument separately so that the range of the function 
is the interval [-7t,+7r] with no singularities.) 


J 


* < 


tan 


-1 


V 


R x 
\ J 


23 2 Aerodynamic Drag 

An aerodynamic drag force results when the moving spacecraft interacts with the ratified 
atmosphere, producing a force that opposes motion. Because this force is, in general, not symmetrical 
with respect to the spacecraft CG, a disturbance torque results as well. (This is discussed in 
section 3.2. 1.3.) The general equation for computing aerodynamic drag is 

F = \ c dPa v2a 


where 

F aerodynamic drag force 
p A atmospheric density 
c D drag coefficient (assumed »2.6) 
v spacecraft velocity 

A projected area of spacecraft normal to spacecraft velocity vector 

The atmospheric drag force is computed for each axis by applying the preceding equation while using 
the velocity (in spacecraft body coordinates) along each axis and by considering the projected area of 
the spacecraft normal to each spacecraft axis. It should be emphasized that this method is an 
approximate one but is sufficiently accurate for the purposes of this model, because the drag torque 
must also be computed, the centroid of this projected area must be specified as well. Therefore, using 
vector notation, 



( F dng}x 

( F dnjy 

(^drapz 


1 

~2 c D?a 


sgn( v s/ C ) JC (Vy c yA x ) 

sgn(v s/c ) > ,(v^ c ) y (Ap 

sgn(v s/c ) z (v^ /c ) z (A z ) 
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where 

drag 

C D 

Pa 


drag force vector (in spacecraft coordinates), lbf 
drag coefficient (dimensionless, equal to 2.6) 

atmospheric density in spacecraft vicinity (from Jacchia’s 1970 model), slug/ft 3 


7 S/C spacecraft velocity vector, ft/sec, transformed into spacecraft coordinates (using attitude 
quaternion) 

A ;C Ay^ l projected area of spacecraft normal to spacecraft x, y, and z axes, respectively, ft^ 

The atmospheric density p A is computed by means of a subroutine that utilizes Jacchia s 1970 model. 
A complete description of this atmospheric density model can be found in reference 2. 


233 Axial Thrust 

2.3.3.1 Thrust vector . — The linear force acting on spacecraft due to axial thrust is desired to be in 
the spacecraft positive x direction; however, because of the thruster-gimbal system and unintentional 
thruster misalignments, there are typically small y and z components of axial thrnst as well. The 
magnitude of the thrust for blowdown systems varies with time and is as given in section 2.33.3. lhe 
thrust magnitude for regulated systems is assumed to be constant at some specified value. The thrust 
vector (in spacecraft coordinates) is expressed as 

F tg cos Y v cos Y z F TF cos cos E z 

f? T = F tg sin Y z + F TF sin e z 

-F tg sin Yy cos y z gimballed ~ F TF sin B y J flxed 

L 7 J thruster thiuster 

where 

F* t resultant thrust vector in spacecraft coordinates, lbf 
F m magnitude of thrust produced by gimballed thruster, lbf 
Fjp total of thrust magnitude produced by fixed thrusters, lbf 
Yy gimbal angle about y axis, rad 

Y z gimbal angle about z axis, rad 

6^ angular misalignment of fixed thrusters about y axis, rad 

angular misalignment of fixed thrusters about z axis, rad 

Positive angles are defined as counterclockwise about the positive axis in question. The fixed thruster 
term uses small angle approximations; the gimballed thruster force is derived in section 23.3.2. 


15 


For COLD— SAT the thrust level of the gimballed thruster Fjq is 0.04 lbf and corresponds to the 
minimum thrust level required to conduct some of the onboard experiments. When an experiment 
requires axial thrust levels higher than 0.04 lbf, F-jj- is commanded to produce a force equal to the 
required thrust less 0.04 lbf. This arrangement provides the required thrust vector control by 
gimballing a single, low-thrust thruster. 


The gimbal angles y y and y z are determined by the attitude control system in a manner to reduce 
attitude error. Refer to section 3.7.2 for a description of this calculation. 


2.33.2 Gimbal model . — The gimballed thruster system consists of a 0.04-lbf hydrazine thruster 
mounted on two single-axis antenna gimbal mechanisms. Gimbal #1 is attached to the spacecraft and 
can rotate as much as ±20° about the spacecraft y axis; gimbal #2 is mounted on gimbal #1 and can 
rotate as much as ±20° about the z axis. The rotation angle of gimbal #1 is designated as y , and the 
rotation angle of gimbal #2 as y z . Therefore, if by defining a coordinate system fixed in the thruster 
itself and assuming that the thruster force vector is directed along the +x axis of the thruster coordinate 
system, the thrust vector produced by the gimballed thruster in spacecraft coordinates can be 
determined by applying this coordinate transformation: 



cos Yy 

0 

sin Yy 

0 

1 

0 

-sin Yy 

0 

COS Yy 


|cos y z 
sin y z 
0 


-sin y z 
cos y z 
0 


0 

0 

1 


0.04 (lbf) 
0 
0 


Multiplying the matrices in the proper order gives the gimballed-thruster thrust vector in terms of the 
gimbal angles, which is the first term of the equation given in section 2.3.3. 1. The desired gimbal 
angles are computed by the rotation model as a function of attitude error, integral of attitude error, and 
attitude rate. 


2.3.33 Thrust decay (blowdown system ).— In order to accurately model the thrust produced by a 
blowdown system, the thrust magnitude Fj = Fjp + Fjq is allowed to vary as a function of hydrazine 
usage in the model by means of a table-lookup algorithm. For a blowdown system F T is a linear 
function of pressure, which drops as hydrazine is expended. The tank is repressurized with helium as 
required during the spacecraft life in order to maintain a workable hydrazine pressure. A typical graph 
of system pressure P ^ as a function of hydrazine expended is shown in figure 2.3. The thrust level is 
related to the system pressure by 

f t ‘ (Ft) 

Miom 


where (Fj)^ is the nominal thrust level, is the system pressure (as given in fig. 2.3), and P nom 
is the nominal system pressure (assumed to be 150 psia for COLD-SAT). The rate of change of 
hydrazine mass (denoted m^,^) is a function of the actual thrust F T as follows: 


rfm N 2 H 4 

Jt 



where I g p is the specific impulse, a function of the system design. If a pressure- regulated hydrazine 
supply system is used, then die effect of system pressure can be eliminated by setting the system 
pressure variable P (as a function of hydrazine mass) to a constant value. 
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2.3.4 Translational Force Resulting From Unpaired RCS Thrusters 

The baseline COLD-SAT RCS system is an uncoupled-thruster system with eight thrusters, which 
produce torque couples only about the jc axis (see section 3.2.3.2). When control thrusters are fired to 
produce torques about the y and z axes with such a system, a linear acceleration will be produced in 
the +jc direction. The COLD-SAT dynamic model accounts for this linear acceleration in the transla- 

tion model as F RCS . 


2.4 Approximation of Spacecraft Mass 

The mass of the COLD-SAT spacecraft varies over the life of the mission because of the various 
experimental and propulsion system fluids used. For the purposes of the model the total spacecraft 
mass is considered to be the sum of three components: 

m S/C ■ m diy + m H 2 + m N 2 H 4 

where is the mass of the spacecraft minus the mass of all expendable fluids, is the mass 

of hydrazine used for all axial thrusters, which is computed as indicated in section 2.3.3 3, and mu is 
the mass of all other expendable fluids onboard, primarily liquid H 2 , which is entered into the model 
as a table of mass versus time. 


2.5 Computation of Spacecraft Position and Velocity Vectors 

The resulting acceleration due to all external forces acting on the spacecraft in body coordinates is 
given as follows: 


(*Ws/C - 


* ^RCS * ^clrag 
m S/C 


This acceleration vector is then transformed into TOD coordinates as follows: 


0 


0 

(^WtOD 

= 9* 

(“cxt^S/C . 


where 

‘ q actual spacecraft attitude quaternion (When the translation model is used alone, the desired 

attitude quaternion is used in place of the actual attitude quaternion because the latter must 
be computed by the rotation model.) 

conjugate of spacecraft attitude quaternion 
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( a ex ,)roD spacecraft linear acceleration due external forces in TOD coordinates 
* signifies quaternion multiplication 

The notation 


0 


represents a quaternion whose scalar part is zero and whose vector part is equal to the vector ct. Refer 
to appendix B for a description of quaternion algebra. 

Thus, the total linear acceleration acting on the spacecraft due to all modeled forces is expressed as 

a*= T + (^xth'OD 

where 


a total spacecraft linear acceleration in TOD coordinates 
g* spacecraft linear acceleration due to Earth gravity in TOD coordinates 
The spacecraft velocity and position vectors can now be computed as 


vtO = [/ afarfr + 

= f vtr)dT + % 

Jt o 


where 


R spacecraft position vector in TOD coordinates 
v spacecraft velocity vector in TOD coordinates 
Rq initial position vector (see section 2.2) 
vj) initial velocity vector (see section 2.2) 
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2.6 Computation of Desired Attitude Quaternion 


2.6.1 Overview 

The COLD-SAT attitude model computes a quaternion that represents the desired attitude of the 
spacecraft The actual spacecraft attitude quaternion is computed by the rotation model. The purpose 
of the desired attitude quaternion is to provide a reference attitude so that the attitu e error (a measure 
of the difference between this desired attitude and the actual attitude) can be computed. The attitude 
control system then uses this information to correct the actual spacecraft attitude so that it is as close 
to the desired attitude as possible. The desired spacecraft attitude options are defined by the alignment 
of certain spacecraft axes with certain base unit vectors, which are defined m TOD coordinates. 

The COLD-SAT dynamic model allows an attitude to be selected from among nine basic options. 
Within each option several suboptions are permitted. The attitude may then be further mo 1 i y 
specifying a fixed or time-varying rotation about any one of the spacecraft axes. Figures . to 
are pictorial representations of each of the attitude options. The computation of the desired spacecraft 
attitude proceeds as follows depending upon the user-specified attitude option, axis alignment, and 
optional spacecraft rotation: 

(1) Computation of the position of the Sun in TOD coordinates based upon the Julian date at the 
start of the simulation and the simulation time 

(2) Computation of two orthogonal unit vectors, which are called base vectors in the discussion that 
follows and are expressed in TOD coordinates. Nine combinations of base vectors, called attitude 
options, are predefined in the model. The desired attitude option is specified by the user at the start of 

the simulation. 

(3) Determination of the spacecraft axes that are to be aligned with each base vector, either positive 
or negative spacecraft axes may be specified. This results in 12 possible attitude venations for each 
pair of base vectors selected. The desired spacecraft axis alignment is specified by the user at the start 

of the simulation. 

(4) Computation of the rotation axis and the rotation angle for the optional spacecraft rotation. The 
rotation axis may be any positive or negative spacecraft axis; the rotation angle is expressed as 


e(r) = it + e 0 

The rotation axis, as well as the values of i and e 0 , is specified by the user at the start of the 
simulation. 

(5) Computation of the desired attitude matrix (direction-cosine matrix) based on steps (1) 
through (4) 

(6) Extraction of the desired attitude quaternion from the desired attitude matrix 
Each of these steps is described in detail in sections 2.6.3 to 2.6.8. 
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2.6.2 Sun Model 


The Sun model computes the position of the Sun in TOD coordinates on a given Julian date, as 
well as a unit vector normal to the plane of the ecliptic, also expressed in TOD coordinates (ref. 3) 
The mean solar longitude with respect to die TOD vernal equinox ^ ns , in degrees, is computed in 
Earth-centered ecliptic coordinates (see appendix C for a description of this coordinate system) as 
follows: 


Km = 280.46° + 0.9856474°(JD 2000 ) 


where (JD 2000 ) is the Julian date referenced to the year 2000. The Sun’s mean anomaly M, in degrees, 
is now computed as 


M = 357.528° + 0.9856003°(JD 2000 ) 

The ecliptic longitude of the Sun in degrees, is computed as 

Sed = + 1.915° sin(Af) + 0.020° sin(2 M) 


The unit vector in the direction of the Sun Ug un and a vector normal to the plane of the ecliptic 
are expressed as follows in ecliptic coordinates: 



c«K5cci) 

singed) 


0 



0 

0 

I 


The same vectors are expressed in TOD coordinates as follows: 



(USuA 


( u SunW) 

U Sun 

( u Sun)y 

= 

(”SunW°s{23.439° -4 x10- 7 (JD 200 o)} 


( u Sun)* 


,(”SunWin(23.439° - 4 x 10- 7 (JD 2000 )j 



( v ecl^x 


0 

II 

& 

( u ecl)y 

= 

-sin(23.439° - 4 x1O" 7 JD 2000 ) 


<”ecl)?_ 


cos(23.439° - 4xlO" 7 JD 2000 ) 
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2.6.3 Base Vectors 


Two orthonormal base vectors, designated v\ and u^, are used to define each of the nine options 
built into the COLD-SAT dynamic model. These base vectors are selected from the normalized 
position, velocity, Sun, and ecliptic-normal vectors or from some mathematical combination of two or 
more of these vectors. Each of the nine attitude options is defined in table 2.1 and illustrated in 
figures 2.4 to 2.12. 

The two base vectors are chosen from the following list of unit vectors as defined in table 2.1, 
depending on which of the nine attitude options is selected: 


u^ un unit vector in direction of Sun (This vector is computed by the Sun model.) 

0^1 unit vector normal to plane of ecliptic (This vector is also computed by the Sun model.) 

ujj unit spacecraft geocentric position vector (This vector is the normalized spacecraft position 
vector R .) 

\? v unit spacecraft velocity vector (This vector is the normalized spacecraft velocity vector v. ) 

4 x K 


v* N unit vector normal to spacecraft orbit plane, | ^ ^ y 
u* unit vector in direction of projection of Sun vector into spacecraft orbit plane, 


% X u Sun 


l«JV x u Sunl 


V N 


o^, unit vector in plane formed by Sun vector and spacecraft velocity vector and normal to 


spacecraft velocity vector, 


X ”Sun 


l»v X U 


Sunl 


x U, 


v*p unit vector in orbit plane and normal to Sun vector. 


% X U Sun 


\VN X u Sunl 


2.6.4 Spacecraft Axis Alignment With Base Vectors 

Once the base vectors are defined, two of the spacecraft coordinate axes are aligned with the two 
base vectors to establish the desired spacecraft attitude. This alignment of spacecraft axes with base 
vectors is controlled by the manner in which the attitude matrix is constructed, as described next. 


21 



2.6.5 Computation of Desired Attitude (Direction-Cosine) Matrix 

The attitude matrix (sometimes called a direction-cosine matrix) defines the transformation of a 
vector from the TOD coordinate system to the spacecraft body coordinate system. The attitude matrix 
is a 3*3 matrix with the property that each column (and each row as well) is a unit vector. If v^ OD is 
a vector defined in TOD coordinates and is the same vector expressed in spacecraft body 
coordinates, the relationship between the two representations is 

vi/C ~ [TOD — * S/C]Vjqd 

where [TOD — ► S/C] is defined as the attitude matrix. 

The attitude matrix is formed by substituting the two base vectors for the selected attitude option 
into the two columns of the matrix that correspond to the two spacecraft axes which are to be aligned 
with these base vectors. The first column of the matrix corresponds to the x axis, the second column 
to the y axis, and the third column to the z axis. If the negative of a given axis is to be aligned with a 
base vector, that vector is negated before substitution into the matrix. Because the base vectors are 
unit vectors by definition, the remaining column of the matrix is found by taking the crossproduct of 
the two vectors in the proper order to form a right-handed system. 


In equation form, in order to align the spacecraft -x axis with the base vector and the z axis 
with the base vector iJj , the attitude matrix would be as follows: 

[-«2 I * ( _ C 2 ) I W |] 

As an example, for the attitude selected for COLD-SAT (attitude option 4, with the spacecraft -x axis 
aligned with and the spacecraft z axis aligned with v^), the attitude matrix would be as follows: 



»N * (-«!*) I 


2.6.6 180° Spacecraft Rotation 

For attitude options 4, 7, and 8, the spacecraft is required to execute a 180° rotation about a given 
axis whenever the Sun crosses the spacecraft orbit plane. This maneuver is necessary to allow the side 
of the spacecraft along which the cryogenic lines are run to remain shaded at all times and to 
minimize the angle between the solar panel normal and the Sun vector. This rotation axis, expressed 
in TOD coordinates, is as shown in table 2.2. The corresponding axis in spacecraft coordinates 
depends upon which axis is aligned with the TOD vector. 

It is undesirable to allow the desired spacecraft attitude to change instantaneously by 1 80° because 
this would result in a sudden 1 80° attitude error. The subsequent response by the attitude control 
system (ACS) to compensate would result in poor transient behavior, including large attitude rates, 
with significant overshoot of the desired attitude and subsequent oscillation about the desired attitude 
with a very long settling time. In order to avoid this problem, the change in desired attitude is 
executed over a finite time tj. 
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The minimum achievable value for tf is 


('/> min = 

where 

/ spacecraft moment of inertia about rotation axis 
t maximum control torque that ACS can produce about rotation axis 

This equation assumes that the spacecraft undergoes a constant angular acceleration at the maximum 
achievable rate during the first half of the 180° rotation maneuver and then decelerates at an equal rate 
for the remainder of the maneuver. 

For any time t f >. the desired attitude matrix is computed as follows: 

A(f) = R 1 (’P)A(r 0/ ) 


where 

A(r) desired attitude matrix at time t, where t is with respect to the start of the rotation (fy) 

A(f(y) desired attitude matrix at instant before time of desired attitude change 

Rjf'F) Euler matrix representing single rotation about appropriate spacecraft axis through angle T 

The an gle T is expressed in radians as 

tp = _Lar 2 for 0 <. t £ -Ify 

2 ^ 

W = 7 T - 0.5 a(t f - t ) 2 for ±t f <t<. t f 

and 

7T 

G 2 

V 

This formulation for T and a assumes a constant angular acceleration for the first 90° of station and 
a constant deceleration for the remaining 90° of rotation. Other time functions for V could be used 
but are not included in the current version of the model. The attitude quaternion is then extracted 
from the desired attitude matrix in the usual manner. 
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2.6.7 Optional Spacecraft Rotation 

For certain attitude options it may be desired to specify a fixed angular offset from the nominal 
attitude or a constant-rate rotation about a given axis (perhaps to keep solar arrays pointing toward the 
Sun for continuously rotating attitude options). The rotation angle is specified as 

e(r) = it + e 0 


where 

t rate of continuous spacecraft rotation, rad/sec 
t simulation time, sec 

e Q constant spacecraft angular displacement, rad 

The desired attitude matrix is formed as follows: 

AM - R ['M]A o0ol 


where 

R(e) Euler rotation matrix of angle e about specified axis 

nominal spacecraft desired attitude matrix (before rotation) 


2.6.8 Quaternion Extraction From Desired Attitude Matrix 


A vector transformation (and thus the spacecraft attitude) can be equivalently represented by means 
of a unit quaternion. Because quaternions are used throughout the model, it is required that the 
attitude matrix representing the desired spacecraft attitude be converted into quaternion form. This is a 
fairly straightforward procedure and is described here. If A is the desired attitude matrix and p is the 
equivalent desired attitude unit quaternion, reference 4 gives A in terms of p as 


2 (p <0> ) 2 + 2 (p <0 ) 2 - 1 


A = 


2 \p%®> + 2/><V 3 > 
2 p<y 3 > - 2 P <v 2 > 


2 p<V 2 > - 2 p%<» 

2 ( p <°>) 2 ♦ 2 ( p < 2 >) 2 - 1 
2 p<»p<» + 2p<°>p< 1 > 


2/»< 1 V< 3 > + 2p%® 

2 pQ>pV> - 2 p<°>p< l > 

2 (/> <0> ) 2 + 2 ( p <3> ) 2 - 1 
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and 


P = 


P 

P 

P 

P 


< 0 >' 

< 1 > 

< 2 > 

<3> 


where is the scalar part of the quaternion. 

Thus, by rearranging terms in the A matrix, the desired attitude quaternion elements can be solved 
for in terms of the attitude matrix elements. 

T = 7r[A] = A n + A 22 + A 33 


Because (p <0> ) 2 + (p <!> ) 2 + (P®) 2 + (p <3> ) 2 = 1 for a unit quaternion, 

<°> = 1 

F 2 


The magnitudes of the three vector elements of the quaternion can be computed as follows: 


I P 


<i> 


+ 


A u M 1 - 7KA) 
2 4 


for 


1,2,3 


The signs of the three vector elements of the quaternion can be computed from the off-diagonal 
elements. The indices i, j, and k are determined in cyclic order so that 

\p®\ * |p 0> |, \p®\ * Ip^l 


Now the vector elements of the quaternion can be determined as 

p® = »g < A kj “ A jk)\P^\ 

p^ = sgn {p^iAji + A^)] |p ^ | 

p^ = sgn[p^(A fa + A^)]|p^| 

where sgn(jc) is the sign function (i.e., sgn(x) returns +1 if the argument x is positive, 0 if the 
argument is zero, and - 1 if the argument is negative). 
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Thus, by using the preceding equations, the four elements of the quaternion can be computed. The 
resulting quaternion is then normalized by means of 



where 

I p\ = 'l(P (0) ) 2 + (/> <0 ) 2 + (P^) 2 + (/> <3> ) 2 


2.7 Miscellaneous Calculations in Translation Model 

2.7.1 Greenwich Hour Angle 

The Greenwich hour angle is the angle between the Greenwich meridian and the TOD X axis. It is 
computed by the following: 


GHA = GHAo + 


where 

GHAo Greenwich hour angle at beginning of simulation, determined by given date and GMT at start 
of simulation (ref. 3) 

^Eaith Earth rotation rate, 4.178074216* 10" 3 deg/sec 

The modulo function is applied to the GHA so that it will remain in the interval (0°, 360°). 


2.72 Axial-Thrust Misalignment Torque 

This torque, which results from the thrust-vector offset from the CG, is computed in the translation 
model. However, because this torque affects the rotational behavior of the spacecraft, it is discussed in 
section 3.2.3.3. 
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BASE VECTOR DEFINITIONS FOR PROPOSED COLD-SAT Al l il UDE OPTIONS 


TABLE 2.1 


Attitude 

option 

Base 

vector, 

tJ; 

Base 

vector, 

Description 

1 

^Sun 

*ad 

First axis is in direction of Sun, and second is normal to ecliptic; spacecraft rotates 
once per year. 

2 



First axis points away from center of Earth along position vector, and second is normal 

3 



Similar to (2), except that first axis points along velocity vector and second is normal 

4 

tJ w or -U w * 


One axis points toward projection of Sun vector into orbit plane, and other is normal or 
negative-normal to orbit plane, in direction opposite the Sun, requiring a 180 rotation 
when Sun crosses orbit plane; spacecraft routes once per year; and 180° rotation occurs 
several times per year for low-Earth orbit. 

5 



Same as (4), except second axis is always positive-normal to orbit plane (i.e., no 180 

6 



One spacecraft axis is in plane defined by Sun and velocity vectors; second axis is 
aligned with velocity vector. 

7 


or-V 

One spacecraft axis is aligned with position vector, other axis is normal to orbit plane; 
spacecraft routes one revolution per ort.it. Includes 180° roUtion when Sun crosses 
orbit plane. 

8 


tf w or 

One spacecraft axis is aligned with velocity vector, a second is aligned with orbit 
normal that is on same side of ort.it plane as Sun. Includes 180° roUtion when Sun 
crosses orbit plane. 

9 

^Sun 


One spacecraft axis is in orbit plane; other points directly at Sun. 


Transition from to occurs when the Sun crosses the spacecraft orbit plane. This maneuver is a spacecraft rotation about a 
given axis and is done to minimize solar panel pointing loss. 


TABLE 2.2—180° SPACECRAFT 
ROTATION AXIS 


Attitude 

Rotation axis vector 

option 

(TOD coordinates) 

4 


7 


8 
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Figure 2.1 .—Spherical TOD coordinate system used in gravity model. 
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Figure 2.2.— Orbital elements. 
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Figure 2.3.— Typical hydrazine pressure profile. Figure 2.4.-Attitude option 1 



Figure 2.5.— Attitude option 2. 
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Figure 2.7.— Attitude option 4. {As Sun crosses the orbit plane, 
the spacecraft Is rotated 180° about the spacecraft axis that is 
in the orbit plane and pointing in the direction of the Sun.) 
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Figure 2.10.— Attitude option 7. (When the Sun crosses the orbit plane, the spacecraft 
rotates 100° about the velocity vector.) 



Sun behind 
orbit plane 


Figure 2.11 . Attitude option 8. (When the Sun crosses the orbit plane, the spacecraft 

rotates 180° about the velocity vector.) 
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3.0 Rotation Model 

The rotation model simulates the attitudinal motion of the spacecraft about its center of gravity. 

This motion is a consequence of various disturbance torques acting on the spacecraft as well as torques 
produced by the attitude control system to counteract the undesirable effects of these disturbance 
torques. The torques that act upon the spacecraft may be categorized into three groups. 

(1) External disturbance torques, which include the gravity-gradient, magnetic, solar pressure, and 
aerodynamic drag torques. These torques are, in general, functions of spacecraft attitude, orbit 
position, velocity, and spacecraft geometry and mass properties. 

(2) Self-induced disturbance torques, which are generated by the spacecraft itself. These torques 
include those produced by thruster misalignments and slosh. They are not directly functions of 
attitude, orbit position, or velocity but are related to spacecraft body rates, attitude error, and spacecraft 
design tolerances. 

(3) Control torques, which are intentionally applied to the spacecraft under the command of the 
attitude control system. These torques are functions of attitude error and attitude rate. In practice 
many types of systems are used to affect attitude control; the COLD— SAT dynamic model includes 

(a) Coupled and uncoupled RCS thrusters 

(b) A gimballed thruster, which provides both axial thrust and attitude control about two 
spacecraft axes 

(c) A reaction wheel, which is used to provide control about a single spacecraft axis 

The sum of all torques acting on the spacecraft causes the spacecraft to experience angular 
acceleration, the integral of which is the angular (body) rate. The angular rate in turn produces 
changes in the spacecraft attitude. In the model the attitude is represented in quaternion form. This 
attitude quaternion is compared with the desired attitude quaternion computed by the translation model 
to arrive at the three-axis attitude error. The attitude error, along with the attitude rate, is used by the 
attitude control system algorithms to produce the appropriate control torques for stabilizing the 
spacecraft attitude. 

This section describes in detail each of the following: 

(1) Simulation of all torques included in the model except for the slosh torque. The slosh torque 
computation, because of its complexity, forms a separate part of the COLD-SAT dynamic model (the 
slosh model) and is discussed in section 4.0. 

(2) Computation of the angular acceleration and the angular rates 

(3) Computation of the actual spacecraft attitude quaternion 

(4) Computation of the attitude error 

(5) Description of the attitude control system algorithms 
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The rotation model also includes the computation of the microgravity accelerations at arbitrarily 
selectable points within the spacecraft. A detailed description of these calculations is deferred to 
section 5.0 of this report. A block diagram of the rotation model is shown in figure 1.3. 


3.1 Interfaces of Rotation Model 

Interfaces between the rotation model and the translation model were given in section 2. 1 . The 
interfaces between the rotation model and the slosh model, as well as input and output quantities of 
the rotation model, are given here. The following quantities are computed by the rotation model and 
passed to the slosh model: 

(1) The actual spacecraft attitude quaternion 

(2) The spacecraft attitude-rate vector 

The slosh torque is computed by the slosh model and passed to the rotation model to compute the total 
disturbance torque. 

The following inputs to the rotation model are user specified at the start of the simulation: 

(1) Spacecraft rotational mass properties, which include moments and products of inertia 

(2) Spacecraft magnetic dipole moment, which is used in computing magnetic torque 

(3) Solar and geomagnetic parameters, which are used in computing atmospheric density and hence 
atmospheric drag torque 

(4) Number, location, and thrust of RCS thrusters 

(5) Reaction wheel moment of inertia 

(6) Initial attitude rates 

In addition, numerous program flag variables can be set to disable portions of the model, including 
attitude control system devices not desired to be simulated. Rotation model outputs typically plotted 
include, but are by no means limited to, attitude errors and rates, microgravity accelerations, gimbal 
angles, and total disturbance torques. 


3.2 Simulation of Torques Acting on Spacecraft 

3.2.1 External Disturbance Torques 

3.2.1.1 Gravity-gradient torque model . — Gravitational and centripetal forces, which maintain the 
spacecraft in its orbit, are equal and opposite only at the spacecraft CG point. At all other points on 
the spacecraft body there is a slight imbalance between these two forces, which results in a net torque 
on the spacecraft. This torque is called the gravity-gradient torque and is a function of the orbit alti- 
tude as well as the rotational mass properties of the spacecraft. 
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Gravity-gradient torque in the COLD-SAT dynamic model is computed by using a spherical Earth 
gravity model as follows: 


t GG “ 


3M 


'S/C 



where 

Tq G gravity-gradient torque acting on spacecraft, ft-lbf 

I spacecraft rotational mass-properties matrix (refer to section 3.3.1), ft-lbf-sec 
l*Eaith gravitational constant, 1 .4076469* 10 16 frVsec 2 

R magnitude of spacecraft position vector, ft 

r spacecraft position vector transformed into spacecraft body coordinates, ft, which are 
transformed from TOD coordinates as follows: 


0 



II 

rs 

* 

'o' 

- R S/C- 


jr. 


*<7 


A typical plot of gravity gradient torque versus time for the COLD-SAT spacecraft over approxi- 
mately one orbit, showing x, y, and z components, is given in figure 3.1. The gravity-gradient torque 
can be seen to exist predominantly about the z axis for this attitude and is of sinusoidal form with a 
period of one-half the orbit period. This plot, as well as all others to follow, are for zero axial thmst 
and without the reaction wheel active (because in the final COLD-SAT design the reaction wheel was 
eliminated). All simulation results in this document are from the same simulation run to allow direct 
comparison between different plots. The attitude option is that selected for the final design. 

3.2.1.2 Magnetic torque. — The spacecraft is assumed to have a net magnetic dipole moment A? s/C 
(in ampere-square meters) that is randomly aligned with respect to the spacecraft coordinate axes. The 
interaction of this magnetic dipole moment with the Earth’s magnetic field produces a torque on the 
spacecraft, which is given as follows: 


Cg = (0.737 17 )[a? s/c x (B^s/c] 

where (®Eaith)s/C is the Earth’s magnetic flux density vector (in tesla) transformed into spacecraft 
coordinates. The conversion factor of 0.73717 is required to convert the computed torque (in newton- 
meters) to English units (foot-pounds-force). COLD-SAT is assumed to be a class 3 spacecraft, with 
a residual magnetic dipole moment of 30 A-m^. 


In order to compute B Elu1h , it is necessary to have a reasonably accurate model of the Earth’s mag- 
netic field. In the COLD-SAT dynamic model the Earth’s magnetic field is modeled as a tilted dipole 

with strength = 8. DIO 22 A-m 2 . The flux density produced by this dipole (B EuA ) M , in 

spherical magnetic coordinates, is given by the following expression: 


37 


( fl Eaith)M “ 


*K 

MO ^Eaith 

2 sin <|> w 

.V 

4tt /?3 

-cos <J> W 


where 

H 0 permeability constant, 47T* 10~ 7 H/m 
R magnitude of spacecraft position vector, m 

<J> W spacecraft magnetic latitude measured from the Earth’s magnetic equator, which is computed as 

4 >m = tan 

where 



R m = [TOD -M]fl 

and R m is the spacecraft position vector expressed in Earth-centered magnetic coordinates (refer to 
appendix C). The matrix [TOD — ► A/] is the transformation from TOD to magnetic coordinates and is 
derived in appendix C. 


The vector as computed previously, is in spherical magnetic coordinates with the R 

component directed away from the center of the Earth and the positive <J> W component directed 

perpendicular to R and toward magnetic north. This vector is converted into the rectangular 
magnetic coordinates as follows: 




B x 

< 

B y 

= 1 


(Earth) w 


sin 


d M cos <J> M 

e A# c ° 8 4 >m 
sin 


■— cos 0 ^ sin <J) ^ 
-sin Q m sin <|> M 
cos <j> w 


h 

B* 




where is the spacecraft magnetic longitude, which is computed as 




tan 


-1 
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The vector flEaith * s now transformed to TOD coordinates: 


(^Eaith^TOD - l TOD M] (^Ewth^M 


Finally, the flux density is transformed into spacecraft coordinates by using the usual quaternion 
transformation as follows: 


0 


k fl Eaith)s/C 


= (q c )* 


(^Eaith^TOD 


*9 


where if is the spacecraft attitude quaternion. 

A typical plot of magnetic torque versus time for the COLD-SAT spacecraft over approximately 
one orbit, showing x, y, and z components, is given in figure 3.2. There is no magnetic torque about 
the x axis because the spacecraft magnetic dipole moment chosen was aligned with the spacecraft 

jc axis. 

3.2.1.3 Aerodynamic drag torque— Aerodynamic drag torque is produced by the effect of 
unbalanced (with respect to the CG) drag forces acting on the spacecraft (refer to section 2.3.2 tor 
computation of the aerodynamic drag force). The following quantities are the coordinates of the 
centroid of the projected spacecraft area along each principal coordinate axis: 

x z ,y z x and y coordinates, respectively, of centroid of projected area normal to spacecraft z axis 

y^x y and z coordinates, respectively, of centroid of projected area normal to spacecraft x axis 

x and z coordinates, respectively, of centroid of projected area normal to spacecraft y axis 

The aerodynamic torque can then be computed as follows: 



T * 


£y(^drag);y ” 

T drag 

T y 

= 




drag 

>x^drag)x ' x y^d ragV 


where Fjng. the aerodynamic drag force, was computed in section 2.3.2. 

A typical plot of aerodynamic drag versus time for the COLD-SAT spacecraft over approximately 
one orbit, showing x, y, and z components, is given in figure 3.3. The apparent ripple in the curves 
results from the fact that the atmospheric density is computed at approximately 1-min intervals to 
reduce execution time. This computation involves iterative routines, and if computed continuously, 
would consume excessive central processing unit (CPU) time. Furthermore, this torque is very small 
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relative to the dominant gravity-gradient torque so that highly accurate computation is not required. 
There is no drag torque about the x axis because the spacecraft as modeled is completely symmetrical 
with respect to the x axis. (CG shift from the nominal CG location is neglected for this computation.) 

3.2.1.4 Solar pressure torque model . — The Sun shining on a surface is equivalent to the surface 
being bombarded by photons; it produces a pressure on the surface. Because, in general, the surface 
area of the COLD— SAT spacecraft upon which the Sun is shining is not symmetrical with respect to 
the CG, this solar radiation pressure will produce a disturbance torque on the spacecraft. 

Highly accurate modeling of the solar pressure torque produced on a spacecraft requires a finite- 
element approach, which is very costly in terms of CPU time. For the purposes of the COLD-SAT 
spacecraft, a much simpler, less accurate approach will suffice, because the solar pressure torque is 
very small relative to the dominant gravity-gradient torque. 

The solar pressure torque modeling approach chosen for the COLD-SAT dynamic model will pro- 
duce correct results only for the attitude option and axis alignment chosen for the COLD-SAT design, 

which is attitude option 4 with the spacecraft z axis aligned with uj , the orbit plane normal, and the 
-x axis aligned with u^, in the orbit plane pointed toward the Sun (refer to section 2.6). 

If the COLD-SAT spacecraft is considered to consist of a cylinder and a truncated cone, as shown 
in figure 3.4, with the CG located at the center of the cylinder, it is apparent that the net torque 
produced by sunlight striking the cylinder and the solar panels (which are perpendicular to the page in 
figure 3.4) will be zero, because both are completely symmetrical with respect to the CG. Solar pressure 
torque will be produced only when sunlight strikes the conical portion of the spacecraft. Since the 
Sun always lies in the x-z plane and the angle p varies between 0° and 41°, the solar pressure torque 
will be zero when the angle 0 is less than 15° because the conical portion of the spacecraft will be in 
shadow. If 0 is greater than 15°, sunlight will strike the conical portion of the spacecraft and a 
negative solar pressure torque will be produced about the spacecraft y axis. 

The value of this torque is computed as follows: As given in reference 5 the solar pressure force 
produced by sunlight striking an object of small surface area dA is given by 


^solar ^solai^ + P)( u Sun " u /») 

where 

Pgoto solar constant, 9.4* 10~ 8 lbf/ft 2 at 1 AU 

p surface reflectivity (p = 1 for highly reflective surface, p = 0 for black surface) 
unit Sun vector 

u; unit vector normal to area dA 

The planar projection of the conical area of the spacecraft in any direction normal to the surface is 
75.44 fr as shown in figure 3.4. The torque moment arm (the distance between the centroid of the 
projected conical area and the spacecraft CG) i is 1 1.3 ft as determined from figure 3.4. Assuming a 

worst-case reflectivity of 1 and realizing that u^ un • = cos 6 = sin(p - 15°), the solar pressure 

torque can be computed by the following equation: 
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where 


^solar ^y^solar’’’ v 


for P £ 15° 


K. 


solar 


solar torque constant, + p)A«, 1.6M0' 4 ft-lbf for COLI>-SAT spacecraft 

solar visibility flag, 1 for visible Sun and 0 for eclipse 


solar beta angle, cos 


-1 


(”^Sun)s/C 


* V 


deg or rad 


unit vector in direction of spacecraft x axis 
unit vector in direction of spacecraft y axis 
(^SuJs/C unit Sun vector ’ transformed into spacecraft coordinates 

It should be emphasized that this solar pressure model is highly specific to COLD-SAT and is 
generally not applicable to other spacecraft. Furthermore, it is only valid for small attitude errors. 

Figure 3.5 is a typical plot of solar pressure torque versus time for the COLD-SAT spacecraft over 
approximately one orbit, showing x, y, and z components. As previously discussed, this torque is 
modeled about the y axis only. In the regions designated A on the plot, the Sun is visible and thus the 
solar pressure torque is nonzero; in region B the Sun is blocked by the Earth and the solar pressure 

torque is zero. 

In order to accurately compute solar pressure torque, it is necessary to determine when the Sun is 
blocked by the Earth. The solar blockage routine in the COLD-SAT dynamic model determines 
whether the spacecraft is in sunlight, penumbra, or umbra and sets the value of the solar visibility flag 
5 to 1 , 1/2, or 0, respectively. Figure 3.6 represents the plane that contains the Sun, the Earth, and 
the spacecraft, showing the umbra and the penumbra. If an angle 0^^ is defined as the angle 
measured from the subsolar point (i.e., the point of intersection between the Earth-center-to-Sun-center 
line and the spacecraft orbit) to the spacecraft position vector, whether the spacecraft is in umbra, 
penumbra, or sunlight can be determined by the following: 



1 

for 

® range < 6 P 

(sunlight) 

S v = 

1 

2 

for 

Q P * ® range < Q U 

(penumbra) 

s v = 

0 

for 

® range > 

(umbra) 


Because the umbra and the penumbra are symmetrical with respect to the Earth-center-to-Sun-center 
line, all angles measured from the subsolar point can be considered to be in the first or second 
quadrant only (i.e., 0° S 0 ^ 180°). Thus, angles measured from either side of the subsolar point are 
considered to be positive. The angle 0 ran g C is computed by 
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® range = * ^Sun) 


where v* R and t^ un are unit vectors as defined in section 2.6.3. The angles Q p and d v are measured 
from the subsolar point to the points where the penumbra and umbra each intersect the spacecraft 
orbit, respectively. They are computed as follows: 

= 90° + 0^ + + 0«it 

where 0^ is the correction for the difference between die Earth’s radius and the Sun’s radius for the 
umbra and is computed as follows: 


and 


®esdu tan 


*Sun _ ^Earth 


°Eaith-: 


Sun 


= 0.26412° 


R Sun solar radius, 2.28346* 10 9 ft 

Earth equatorial radius, 20 925 656 ft 
f>Eaith-Sun Earth-Sun distance, 1 AU or 4.908066* 10 11 ft 
The correction for atmospheric refraction, 0^ (in degrees), is given as follows: 


0„ = i.i33 - O-seec" 017 ** 11 


where JCjjj is the spacecraft altitude in nautical miles. This term is essentially constant at 1.133 for 
altitudes above 100 nautical miles. The altitude correction 0^ is computed as follows: 


For the penumbra, 


©alt = cos 


lf^Ewth) 



= 90 ° - ©csdp + e ar + e at 


where 0 esdp is the Earth and Sun radii correction term for the penumbra, given as 


0 csdp = ton_1 


^Sun + ^Earth 


■^Earth- 


Sun 


= 0.26901° 


and 0^. and 0 a]t are the same as for the umbra. 
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3 22 Slosh Torque 


The torque acting on the spacecraft is computed by the slosh model, which is discussed in detail in 
section 4.0. 


3.23 Thruster-Generated Torques 

3.2.3.1 General . — Any thruster produces a translational force on the spacecraft in the direction of 
the thrust vector. If the thrust vector does not pass through the spacecraft CG, it will also pro uce a 
torque. A spacecraft thruster is classified as an RCS thruster if its purpose is to produce a control 
torque on the spacecraft; as a fixed axial thruster if its purpose is to produce a translational force on 
the spacecraft; or as a gimballed axial thruster if its purpose is to produce both a translational force 
and a control torque. However, because of thruster misalignments and spacecraft CG shifts, m 8 enera 
all of these classes of thrusters will produce both rotational torques and translational forces, and the 
undesired force or torque will constitute an additional disturbance on the spacecraft translationa or 
rotational dynamics. All three of these thruster classes are present on the COLD-SAT spacecraft and 
are discussed in detail in subsequent sections. 

The purpose of the RCS is to provide control torques on the spacecraft in order to maintain the 
desired spacecraft attitude. The primary advantages of an RCS over other attitude control schemes are 
relatively low cost, less complexity, high reliability, and ability to control both cyclical and secu ar 
disturbances for axial-thrust and no-axial-thrust conditions. The principal disadvantage of such a 
system is that the discontinuous nature of this type of control degrades the experimental microgravity 
environment (see section 5.0). An RCS can be classified as a coupled or uncoupled scheme 
depending on whether or not the torques produced are pure couples. An uncoup ed RCS scheme, 
which is discussed in detail in section 3.2.3.2, will always produce a translational force on the 
spacecraft regardless of whether thruster misalignments are present. This translational force will resu 
in additional microgravity disturbances and will, over time, affect the spacecraft orbit However, this 
uncoupled RCS scheme is advantageous from a cost or complexity point of view, and despite its 

disadvantages, was chosen as the baseline RCS for COLD-SAT. From a simulation point of view, 

however, the two types of systems are the same, except that the thrust produced by certain thrusters is 
set to zero when simulating the uncoupled system. 

A fixed axial thruster is designed to produce a translational force on the spacecraft for the purpose 
of inducing a controlled microgravity acceleration inside the experiment tanks. However, because of 
CG shift and thruster misalignments, a torque will also be produced as previously discussed. This 
torque constitutes a disturbance that must be counteracted by the attitude control system. 

A gimballed axial thruster is a single low-thrust device attached to a bi-axial gimbal mechanism. 
This device serves two purposes: First, it augments the thrust produced by the fixed thrusters m order 
to maintain the desired acceleration of the spacecraft; and second, it constitutes part of the attitu e 
control system because the orientation of the thrust vector produced by this thruster with respect to the 
CG can be controlled, thus providing a control torque on the spacecraft. 

The torque produced by any thruster X* can be computed by the simple crossproduct rule. 


-T= F** F 
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where 


F thrust vector produced by thruster 

r vector from location of thruster (i.e., point of application of thrust vector to spacecraft rigid 
body) to spacecraft CG 


For each thruster the values of F* and F are given in tabular form as part of the user-supplied model 
input data. The preceding equation is used to compute the torque for each thruster that is turned on. 
The torques from each individual thruster are added to obtain the total thruster-produced torque acting 
on the spacecraft. 

3.2.3.2 RCS thruster torques . — The COLD— SAT dynamic model can model any desired number of 
RCS thrusters, whose steady-state thrust level, orientation, and location with respect to the CG can be 
specified by the user. For COLD-SAT two different RCS schemes were considered: 

(1) A coupled RCS system, which utilizes twelve 0.06-lbf thrusters located as shown in figure 3.7. 
The thrusters are always fired in pairs so that the resulting torques are always pure couples (i.e., there 
is no translational force acting on the spacecraft when a thruster is fired). 

(2) An uncoupled RCS system, which utilizes the same eight rear thrusters; however, the four front 
thrusters (2, 4, 6, and 8) are eliminated. This scheme results in considerable cost savings, because the 
four front thrusters and the associated hydrazine lines and line heaters running from the tanks in the 
rear of the spacecraft are eliminated. However, the torques produced about the y and z spacecraft axes 
are no longer couples (i.e., there is a translational force produced in the +x direction whenever a y axis 
or z axis RCS thruster is fired). Also less control torque is available about the y and z axes because 
only a single thruster is fired instead of the two thrusters fired with the coupled scheme. Note that the 
torques produced about the x axis are still pure couples. 

In addition, it was assumed that all RCS thrusters are misaligned with respect to the spacecraft 
coordinate axes by 1°, so that when a control torque is required about a given axis and the appropriate 
thruster (or thrusters) is fired, disturbance torques will be produced about the other two spacecraft axes 
as well. These misalignments are included in the COLD-SAT dynamic model and are chosen such 
that the resulting disturbance torques are maximized. The following data are entered to simulate any 
RCS system in the COLD-SAT dynamic model: 

(1) The number of RCS thrusters 


(2) A table of the steady-state thrust vector produced by each thruster, including any misalignments 

(3) A table of the coordinates of each thruster in the system (see section 2.0). The distance 

of the CG from the plane containing the back thrusters is entered into the axial thrust model. CG 
shift from nominal along the y and z axes is neglected in RCS torque calculations. 

For COLD— SAT the number of thrusters is set to 12 for either the coupled or uncoupled scheme. The 
thrust vectors produced by each thruster are as given in table 3.1. In order to simulate the uncoupled 
scheme, the thrust vectors produced by thrusters 2, 4, 6, and 8 are set to zero. 

The coordinates of each thruster in the x',y',7! system are specified by the user. The coordinates of 
the thrusters relative to the spacecraft CG can then be readily computed by specifying the spacecraft 
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CG location relative to the origin of the x'./.z' system. These torques are computed only once at the 
beginning of the simulation, in order to reduce CPU time; therefore, subsequent movement of the CG 
will not be accounted for in the RCS simulation. However, this is not a serious model limitation 
because simulation of the full-up model rarely, if ever, is run for more than a few hours, over which 
time the CG can be expected to move very little. Because the y and z CG offsets, (y y ) an (z Z), 
are very small, these offsets are ignored in computing RCS torques. These thruster locations (table 3.2) 
assume that the front thruster ring is 286.91 inches from the back ring and that the thruster ring radii 
are 48 inches for the back thruster ring and 20 inches for the front thruster ring. 


From these coordinates a vector (^ cs ), is formed from the thruster to the CG for each thruster as 
follows: 


('RCs)i “ 


“(^ RCsX “ (*CG)i 
-(/rCs)i 
RCS)/ 


where 


(Acs)/’ (/Res),- coordinates of ith RCS thruster in spacecraft xf,y\ Z system as given in 
(*RCS>i toble 3 - 2 

x ' spacecraft CG x coordinate in x 1 , y', z' spacecraft coordinate system (i.e., distance 

of CG from rear thruster ring) 

Recall that / CG and z' CG are assumed to be zero and thus are not included in the computation of 

(^RCsV 


Next, the torque produced by each thruster is found by 

(TrCs)/ = ( S trans)«[('RCs)i x ( F RCs)i] for 1 = 1,12 


where 

( S trans )i factor between 0 and 1 representing startup and shutdown transient of thruster i as 

discussed in section 3.2.3.2 


(^RCs)< thrust vector produced by ith RCS thruster, as given in table 3.1 

On the basis of the thruster locations and thrust vectors given for COLD-SAT, the resulting torque 
and translational force produced by each thruster are computed in the model. The information given 
in tables 3.3 and 3.4 is presented for informational purposes and illustrates the torques and 
translational forces for a typical simulation run, for thrusters that are perfectly aligned with ® 
spacecraft primary axes (table 3.3) and for thrusters with selected 1° misalignments (table 3.4). The 
direction of each thruster misalignment is chosen so that the resulting disturbance torques are 
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maximized. For both tables x' cc is assumed to be 85.5 inches. Thrusters are always fired in pairs 
(1&2, 3&4, 5&6, 7&8, 9&10, 1 1&12) as discussed in later sections. 

The design of RCS thrusters is such that when a thruster is fired, it will remain “on” for a 
predetermined minimum amount of time, designated r min ^ This is done in order to increase the 
efficiency of the RCS. Furthermore, each thruster exhibits startup and shutdown transients, denoted t m 
and trf, respectively. In reality, the startup transient is a second-order underdamped step response and 
the shutdown transient is a decaying exponential, but in order to reduce unnecessary overhead, it is 
reasonable to model these transients linearly as shown in figure 3.8. These transients are modeled by 
multiplying the steady-state RCS torque for thruster i, as computed previously in this section by a 
factor (St,**), of 0 to 1. If thruster i is in a full-on or full-off state, then S (tram) = 1 or S (tnins) =0, 
respectively; if thruster i is in a startup or shutdown transient, then 0 < Sfj lam y < 1 . 


3.23.3 Axial-thruster-induced torque . — The axial thrust vector F T discussed in section 2.3.3 is the 
vector sum of the thrust vectors that are produced by the fixed and gimballed thrusters, Fj- f and Fjq, 
respectively. If this thrust vector does not pass directly through the CG, a torque on the spacecraft 
will be produced. This torque comprises both a disturbance torque and a control torque. If the 
gimballed thruster were not present, or if the gimbal mechanism were to fail, the resultant thrust vector 
would not, in general, pass through the CG, and thus a disturbance torque would result whenever the 
axial thrusters were fired. With the gimballed thruster in operation, however, the direction of the 
gimballed-thruster thrust vector can be used to provide an attitude control torque of sufficient magni- 
tude to counteract all expected disturbance torques, including that produced by the thrust-vector-CG 
offset of the fixed thrusters. This attitude control scheme will eliminate virtually all RCS thruster 
firings during periods of axial thrust and result in an excellent experimental microgravity environment. 

The torque produced by the axial thrusters (fixed and gimballed) is given as 


t t 


r** Fj 


where 

F r thrust vector computed in section 2.3.3 
r * vector from origin of x ! sc , y' sc , z'j.q coordinate system to CG 

The total thrust vector resulting from fixed and gimballed thrusters is assumed to pass through the 
origin of the x's/q, y's/o ^ S/C coordinate system. It is defined as follows: 


r = 


-^CG 

y ' cg 

^CG 


where x'qq, y' CG , and z! CG are the coordinates of the spacecraft CG. 
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The gimballed thruster is assumed to be driven by a stepper motor, which is modeled as a 
second-order system (with quantized input) as follows: 


M(s) = 


to., 


5 2 + 2£(o m s + 


It is assumed that X, = 0.7071 and to„ = 0.00584, which results in 5-percent overshoot and a time to 
first peak (for a step command) of about 12.6 min. The input of M{s) is quantized with a step size of 
0.07815 deg/step. 


3.2.4 Reaction Wheel Torque 

A single-axis reaction wheel has been considered for COLEV-SAT in order to improve the 
microgravity environment during periods of zero axial thrust. The dominant disturbance torque 
encountered under these conditions is the gravity-gradient torque. This torque appears predominantly 
about the spacecraft z axis and is cyclical. Thus, a single, modestly sized reaction wheel would be 
able to remove this disturbance torque over long periods of time without saturating and would 
eliminate the need for nearly all RCS firings during periods of zero axial thrust. A block diagram of 
the reaction wheel model in rotational dynamics is shown as part of figure 1.3. 

The single-axis reaction wheel in the model includes static and viscous friction. Viscous friction is 
computed as follows: 




where 

F v viscous friction coefficient 

reaction wheel angular speed, rad/sec 
The torque developed by the motor is given as 

T dev - *m V A 

where 

K m DC motor torque constant, ft-lbf/V 
V A DC voltage applied to motor, V 
Thus, the total torque is 

T total = T dev - ^Jv 
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A switch state S c ^ is used to indicate if the wheel is stuck because of static friction. If ■V s 
positive, the wheel is spinning in a positive direction; if 5^- is zero, the wheel is not moving; and if S c y 
is negative, the wheel is spinning in a negative direction. The value of 5^ is determined as follows: 


S c/(' + ) = +1 for ^c/ = °) and ( T dev > T c) or ( S cf = +1 ) and ( W rw > TOL ) 
V' *) = 0 for (S# = 0) and (| Tdev | < t c ) or (S# * 0) and (|wj < TOL) 
S c y(t + ) = -1 for (S cf = 0) and (r dcv < -t c ) or (S^= -1) and (io w < -TOL) 


where 

t c static friction torque that is constant based on physical characteristics of reaction wheel, ft-lbf 
TOL very small positive number (TOL « 1), typically on order of 10 
The rate of change of wheel angular speed is given by 

“nv=^ for 5^*0 

/ >V 

^rvw “ 0 for S cf= 0 


where I w is the wheel moment of inertia; thus, 

w nv = X' UrJWM 

The resulting control torque on the spacecraft when the wheel is spinning is 



If the wheel is stuck because of static friction, the control torque on the spacecraft is zero, even though 
T total ma Y nonzero. 

The reaction wheel angular momentum (in inertial coordinates) must be accounted for in the 
spacecraft equations of motion as follows: 
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0 
0 

Av^nv + 

where w z is the z component of the spacecraft angular body rate vector. 



3.3 Spacecraft Rotational Dynamics 

33.1 Rotational Mass Properties 

Spacecraft rotational mass properties are in the form of a moment-of-inertia matrix, defined as 


*MP 


A Ary Acz 
Ary A’ Az 
Az A'Z A 


where /_, L, and I 7 are the moments of inertia, defined as 

x y 


- J > 2 * 'y - /sc * 1 2 * z ^ 


and Ijy, l xz , and l yz are the products of inertia, defined as 


- -L iy dm - '» ~ ~Lc xz *"• 


The matrix I MP is a user-specified input to the model. 


'^L ix2ty2)dm 
- -L yz dm 


333 Computation of Angular Body Rates 

The rotational dynamic equation of motion is a first-order vector differential equation that allows 
one to compute the angular rates as a function of the spacecraft mass properties and applied torques. 
Euler’s equation for general rigid-body dynamic motion is as follows: 



h 


where 

h angular momentum, ft-lbf-sec 

f* applied torque, ft-lbf 

c? spacecraft angular rate vector, rad/sec 
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Because h = I MP to\ where I MP is a constant matrix, Euler’s equation can be expressed as 

Imp*- & * Omp® 

Both sides of the equation can then be premultiplied by (I^p) -1 to obtain the rate equation 


r-l 


w 


Expanding gives 
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The sum of all torques, both disturbance and control, acting on the spacecraft f* is defined as 


= x = r 


disturb control 


where are disturbance torques acting on the spacecraft: 

'disturb ~ T GG + 'drag + 'solar + 'mag + 'slosh 

an< * 'control are torques produced by the attitude control system and the axial thrusters: 

T control = t RCS + T RW + J T 

The angular rates are then found by integrating the dynamic equation of motion as follows: 

tot/) = (r)dj + wj, 
where o5q is the initial angular rate vector. 
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A typical plot of attitude rate versus time for the COLD-SAT spacecraft over approximately one 
orbit, showing x, y, and z components, is given in figure 3.9. This figure is discussed in more detail 
in section 3.7.1. 1 . 


3.4 Spacecraft Attitude Computation 

The spacecraft attitude represents the orientation of the spacecraft axes (x, y, and z) with respect to 
the inertial (TOD) frame (X, Y, and Z). This orientation can be conveniently specified in terms of a 
quaternion, which represents the rotation required to transform the TOD axes to be coincident with the 
spacecraft axes. The equation describing the time evolution of the quaternion is a first-order vector 
differential equation that specifies the rate of change of the attitude quaternion in terms of the 
spacecraft angular rates and the instantaneous spacecraft attitude as follows: 



where is a time-varying, skew-symmetric matrix formed from the spacecraft angular rates as 
follows: 


0 

-CO 

jc 
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-(0 
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(0 
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-CO 
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K 

<0 
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-to 

w x 
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and 


to = 


W v 


to,, 


to. 


are the angular body rates about the spacecraft x, y, and z axes, respectively, and 


<7 = 


9 

<7 

_<7 


<o>- 

<1> 

<2> 

<3> 


is the attitude quaternion. The quantity q is the attitude quaternion rate, which is integrated to obtain 
the attitude quaternion at any time t. 
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q (') = j['?(t)rfT + q Q 


where q 0 is the initial attitude. 


3.5 Attitude Error Computation 

The attitude error is a measure of the difference between the spacecraft’s actual attitude q and the 
desired attitude computed by the attitude model p (see section 2.6). It is expressed in terms of angles 
about the spacecraft x, y, and z axes, which are denoted 0 X , 0 y , and 0 r respectively. The error 
quaternion (which represents the transformation between die desired spacecraft axes and the actual 
spacecraft axes) is the product of the conjugate of the actual spacecraft attitude quaternion and the 
desired attitude quaternion. 


= <* C )'P 


The three-axis attitude error is defined as 


where 


- -1 


r <oi 

e* 



0 y 

= 

E4 ?err 

0 Z 





L e< ?crrJ 


e = — 1^— for <j> * 0 

sin 4> 

e = -2 for <f> = 0 


and 


• -i <o> 

♦ = cos 

Angle <t> must be in the first or fourth quadrant, so that if (<|> > 7t/2) then let <j> = <j> - it. 

A typical plot of attitude error versus time for the COLD-SAT spacecraft over approximately one 
orbit, showing x, y, and z components, is given in figure 3.10. This figure is discussed in more detail 
in section 3.7. 1.1. 
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3.6 Attitude Determination System and Error and Rate Quantization 

For the purposes of this model the attitude determination system is assumed to provide perfect 
knowledge of the actual spacecraft attitude errors and attitude rates. For COLD-SAT the actual 
attitude determination system consists of rate-integrating gyros, which are updated by information from 
redundant Sun sensors and horizon sensors. These systems, and any errors that they may introduce 
into knowledge of the attitude errors and rates, are not included in the COLD-SAT dynamic model. 
However, the effect of digital quantization of both the attitude errors and rates is included in the 
simulation (although all floating-point numbers stored in a digital computer are quantized, the 
quantization discussed in this section is many orders of magnitude “coarser” than that encountered 
with a 64-bit word length and is meant to represent that found on a spacecraft computer). For 
COLD-SAT the following values are used: 

(1) Attitude-error quantization, 1.43* 10 4 deg/step 

(2) Attitude-rate quantization, 1.43*10 4 deg/sec-step 

These values were chosen to be consistent with the numerical accuracy of the spacecraft computer 
selected for COLD-SAT. 


3.7 Attitude Control System Algorithms 

3.7.1 RCS Thruster Control Algorithm 

RCS thrusters can have one of two possible command states: full on or full off. The decision to 
fire a given thruster is determined by comparing, for each of the three axes, the spacecraft attitude 
errors and angular rates against a set of switchlines on the phase plane (6, 0). Depending on the 
results of these comparisons, a table lookup is implemented to determine which thrusters, if any, need 
to be commanded “on” in order to produce the required corrective torques. 


3.7.1.1 Switchline control . — Reaction control systems are commonly controlled by means of 
switchlines For COLD-SAT the switchline configuration used consists of a pair of parallel straight- 
line segments plotted on the state plane of attitude rate versus attitude error that divide the plane into 
three control regions. One of the switchlines intercepts the ordinate (attitude-rate axis) at a positive 
attitude rate; the other intercepts it at a negative attitude rate. Both have a negative slope and 
corresponding segments have equal slopes. One switchline crosses the abscissa (attitude-error axis) at 
the maximum allowable positive attitude error, the other crosses the abscissa at the maximum 
allowable negative attitude error. The slopes and intercepts of the switchlines for a given axis are 
computed from the RCS-produced maximum acceleration and spacecraft moment of inertia about that 

axis. 

A typical set of switchlines, such as those used for COLD-SAT, is shown in figure 3.1 1. The 
switchlines divide the (0, 0) plane into three control regions: a zero-torque region in the vicinity of 
the origin, a negative-torque region above and to the right of the right switchline, and a positive-torque 
region below and to the left of the left switchline. This scheme is implemented for each of the three 
spacecraft axes. 

A typical phase-plane plot of the z-component of attitude rate versus attitude error is given in 
figure 3.12 for approximately one cycle (~65 min). The presence of the switchlines is clearly evident 
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on the right and left sides of this figure. The regions A, B, and C are marked as in figures 3.9 
and 3.10. As shown in figures 3.9 and 3.10 the primary attitudinal activity occurs about the z axis 
because the dominant disturbance torque (in the absence of axial thrust) is the gravity-gradient torque, 
which exists predominantly about the z axis. Recall that this torque is sinusoidal with a period of 
one-half the orbit period. Thus, when the gravity-gradient torque is positive, it produces a positive 
angular acceleration about the z axis, which causes both the attitude error and the attitude rate to 
increase (region A in figs. 3.9, 3.10, and 3.12). This continues until an RCS switchline is reached. 
When the switchline is crossed, the appropriate RCS thrusters are fired, quickly reducing the attitude 
rate to a very small value. The thrusters are fired about five times per minute, with thrust duration of 
0.2 sec (r min on ) per firing, to oppose the gravity-gradient torque (region B in figs. 3.9, 3.10, and 3.12). 
The phase trajectory hovers around the 1° (0 rev/day) point until the gravity-gradient torque becomes 
negative. The spacecraft then experiences an angular acceleration in the opposite direction toward the 
opposite switchline (region C) until that switchline is crossed and the appropriate RCS thrusters are 
fired. This limit-cycle behavior continues indefinitely as long as there is no axial thrust. 

3.7.1.2 Computing switchUnes . — A switchline is specified by means of a slope and an intercept 
with the attitude-rate axis. The slope and the intercept are chosen so that a positive or negative 
attitude rate outside the “zero torque” region will result in a continuous RCS thruster firing that will 
reduce this attitude rate to zero along a parabolic trajectory in the state plane. Consequently, the slope 
and intercept for a given switchline are functions of the moment of inertia and the steady-state RCS 
torque about that axis. 

For a given maximum-attitude-error specification, steady-state RCS torque, and moment of inertia 
for each axis, the switchlines for each control axis are computed as follows: 



Aslope)/ ~ 


«« 


( 9 ,), 


max 


( L imcpt) f = 


max 


where i = x, y, z, and 

T, magnitude of steady-state, RCS-produced torque about ith axis 

/,• moment of inertia about ith axis from mass properties matrix 

a,- resulting steady-state angular acceleration about ith axis 

(0/)m«x magnitude of maximum allowable attitude error about ith axis 

3.7.1.3 Thruster commanding . — In the COLD-SAT model the decision of which thrusters to fire 
proceeds in two steps. First, the attitude error and the attitude rate are compared against the 
switchlines to determine the “required torque region” (i.e., either positive, negative, or zero, as 
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discussed in section 3.7. 1.1). This is done for each spacecraft axis. Next, the model searches through 
a thruster control table to determine which thrusters are to be turned on. The thruster control table is 
created by the user and reflects the physical arrangement of control thrusters. Table 3.6 is the thruster 
control table used for the COLD-SAT spacecraft. For example, if the attitude error and attitude rate 
about the x axis indicate a point on the x switchline plot in the positive-torque region while the 
attitude error and attitude rate about the y axis give rise to a point on the y switchline plot in the 
zero-torque region and the attitude error and attitude rate about the z axis indicate a point on the 
z switchline plot in the positive-torque region, this combination corresponds to the eleventh line in 
table 3.6. Therefore, for this case thrusters 7, 8, 11, and 12 are fired until there is a change of torque 
region on one or more of the three axes. 


3.7 2 Gimballed Thruster Control Algorithm 

A linearized, small-angle, block diagram representation of the gimballed thruster and the spacecraft 
rotational dynamics is shown in figure 3.13. This representation is reasonably accurate for gimbal 
angles as large as about 20° and thus can be utilized for the purposes of designing the control law. It 
should be emphasized, however, that the actual model simulates the complete nonlinear dynamics of 
the gimbal mechanism, including the fact that the two gimbal angles are not completely independent 
(because one gimbal mechanism is mounted on the other). 

A control law was designed to command the gimballed thruster based upon the spacecraft attitude- 
error and attitude-rate vectors. The control law should meet the following goals: 

(1) The control law must result in a stable closed-loop system for all expected spacecraft 
perturbations. 

(2) The control law must result in a reasonably fast transient response, with a time to first peak of 
approximately 10 to 13 min. 

(3) The control law should not result in excessive overshoot or oscillation. 

(4) The control law must result in acceptable operation for all expected thrust levels. 

The following should be noted concerning the system to be controlled: 

(1) The open-loop system has two poles at the origin (two integrations) as a result of the spacecraft 
dynamics. Any proportional control scheme acting upon the attitude error will result in closed-loop 
poles in the right-half s plane, and thus an unstable system will result. 

(2) An integration in the forward loop is required to have zero steady-state error between the 
attitude error and the gimbal angle. The system as it is, with no integration in the forward path, will 
require a nonzero steady-state attitude error in order to drive the gimbal angle to a nonzero value. 
Because a zero steady-state attitude error is desired, the controller must have one integration (i.e., a 
pole at the origin of the s plane). 

(3) This additional integration in the forward path results in three poles at the origin and thus in an 
even more unstable system. 
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The following approach is used in order to solve this problem. 

(1) If an additional attitude-rate feedback loop is added, one of the poles at the origin is moved into 
the left-half s plane. 

(2) An additional controller zero at some point in the left-half plane near the origin will cause the 
closed-loop root locus to circle about this zero, so that the remaining two closed-loop poles will move 
into the left-half s plane for increasing controller gain. A finite zero and integration can be provided 
by a proportional-integral (P-I) control law. 

(3) Thus, the control law should be a P-I scheme acting on the attitude error, along with a constant 
attitude-rate feedback. Root locus analysis of the system indicates that the following is a suitable 
control law: 


C(s) 


= 2.5(s + 0.0037) = 25 + 


0.00925 


= -400 


Note that this control law is not necessarily optimal but will result in acceptable response for all 
expected conditions. 


3.73 Reaction Wheel Control Algorithm 

The control law used for the reaction wheel is a proportional law that acts on the attitude error 
about the z axis and gives the voltage applied to the reaction wheel drive motor 

V A = 

where 

K p proportional control constant, 300 V/rad 

0, z component of spacecraft attitude error (see section 3.5), rad 
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TABLE 3.1-STEADY-STATE THRUST VECTORS FOR 


COLD-SAT RCS THRUSTERS 


Thruster 

Thrust vector, including 1° misalignments, lbf 

x component 

y component 

z component 

1 

+0.06 cos 1 0 

+0.06 sin 1° 

0 


-0.06 cos 1° 



1 


+0.06 cos 1° 





-0.06 cos 1° 



k 


+0.06 cos 1° 

0 

+0.06 sin 1° 

e 

-0.06 cos 1° 



1 


+0.06 cos 1° 





-0.06 cos 1° 



* 


+0.06 sin 1° 

-0.06 cos 1 0 

0 




+0.06 cos 1° 

1 

mSm. 



-0.06 cos 1° 


mm 



+0.06 cos 1° 

l 


“For an uncoupled RCS system, thrusters 2, 4, 6, and 8 are 
deleted. 


TABLE 3.2.— COLD-SAT RCS 


THRUSTER COORDINATES 


Thruster 

Thruster coordinates in 
spacecraft system, 

in. 


X* 

y' 

z‘ 

1 

0 

n 

0 


-48 

2 

+286.91 



+20 

3 

0 




+48 

4 

+286.91 



-20 

5 

0 


+48 

0 

6 

+286.91 

-20 

1 

7 

0 


-48 


8 

+286.91 

+20 

\ 

9 

C 



) 

-48 

10 





+48 

11 





+48 

12 





-48 
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TABLE 3.3 — RCS TORQUES AND FORCES FOR STEADY-STATE COUPLED SYSTEM 


WITH MISALIGNMENT 


Pair 

RCS thruster 

< 

Steady-state torque, 
ia-lbf 

Translational force, 
lbf 

About x axis 

About y axis 

About z axis 

Along x axis 

Along y axis 

Along z axis 

1 

1 


-2.88 



+0.06 




2* 


-1.20 

— 

-.06 



2 

3 

— 

+2.88 

— 

+.06 





4* 


♦ 1.20 

— 

-.06 





3 

5 



-2.88 

+.06 






e 


— 

-1.20 

-.06 




4 

7 

— 

— 

♦2.88 

+.06 






8* 

— 

— 

♦ 1.20 

-.06 





5 

9 

-2.88 

— 

♦5.13 

— 

-0.06 





-2.88 

— 

-5.13 



+.06 


6 


♦2.88 

— 

♦5.13 



-.06 




♦2.88 

— 

-5.13 

— 

+.06 

— 


•For uncoupled RCS system, thrusters 2, 4, 6, and 8 are deleted. 


TABLE 3 4 — STEADY-STATE RCS TORQUES WITH 1° ANGULAR MISALIGNMENTS 


[Boldface type indicates desired torque produced; normal type indicates misalignment torque.] 


Pair 

RCS 

thruster 

Steady-state torque, 
in. -lbf 

Translational force, 
lbf 

About x axis 

About y axis 

About z axis 

Along x axis 

Along y axis 

Along z axis 

1 

1 

+0.0502629 

-2.87956 

-0.0895309 

+0.0599991 

+0.0010471 


) 


2* 

-.0209429 

-1.19982 

+.210895 

-.0599991 

i 



2 

3 

-.0502629 

+2-87956 

-.0895309 

+.0599991 





4“ 

+.0209429 

+1.19962 

+.210895 

-.0599991 

l 



3 

5 

+.0502629 

+.0895309 

-2-87956 

+.0599991 

0 

+0.0010471 


e 

-.0209429 

-.210895 

-1.19962 

-.0599991 

1 



4 

i 

-.0502629 

+.0895309 

+2-87956 

+.0599991 





8* 

+.0209429 

-.210895 

+1.19962 

-.0599991 

t 



5 

9 

-2-87956 

-.0502629 

+5.12922 

+.0010471 

-.0599991 

0 



-2-87956 

+.0502629 

-5.12922 

1 

+.0599991 



6 


+2-87956 

+.0502629 

+5.12922 


-.0599991 





+2.87956 

-.0502629 

-5.12922 

1 

+.0599991 




‘For uncoupled RCS system thrusters 2, 4, 6, and 8 are deleted. 


TABLE 3.5.— RCS THRUSTER STATES 


Thruster 

state 

Description 

Value of S._ 

0 

Thruster off 


1 

Thruster in startup transient 

0 <5 tang < 1. 

^tnau increasing linearly 

2 

Thruster foil on 

5 tn»n. = 1 

-1 

Thruster in shutdown transient 

0<5 tnui# < 1, 

S tmu decreasing linearly 
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Gravity-gradient torque, tn.-lbf 




(b)y axis. 



(c) z axis. 

Roure 3.1 .—Typical plot of gravity-gradient torque versus time 
for COLD-SAT spacecraft. (Orbit period, 105 min.) 
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(a) x axis. 
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(b) y axis. 



(c)z axis. 

Figure 3.2. — Typical magnetic torque plot. (Spacecraft dipole moment 
aligned with x axis.) 
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■5 l— 


(a) x axis. 


2x1 0" 8 



6x1 0 -6 



(c) z axis. 

Figure 3.3.— Typical aerodynamic drag torque. (Apparent ripple is 
due to the fact that atmospheric density is computed at intervals 
of 1 minute to reduce CPU time.) 





Figure 3.4.— Solar pressure torque computation. (The y axis Is perpendicular to and Into the page.) 
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0 20 40 60 80 100 120 


Time, min 
(c) z axis. 

Figure 3.5. — Typical solar pressure torque plot. 

r~ Point of intersection of spacecraft 


Instantaneous ' °rblt and penumbra 



Figure 3.6. — Sun blockage. 
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Back 

thruster 



O Back-facing thruster 
(force in +x direction) 

% Front-facing thruster 
(force in -x direction) 

Front thruster ring 

Figure 3.7.— RCS thruster location for COLD-SAT spacecraft. (For an uncoupled system the front thruster ring is deleted.) 


t rcs 



-f. 


-trr 


*min on 

Startup I Steady-state 

transient | operation 


Shutdown | 
transient 


! 


I 


Thruster 

commanded 

on 


Thruster 

commanded 

off 


Figure 3.8.— Typical RCS thruster pulse, including transients. 
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(b) y axis. 

Region Region Region 



0 20 40 60 80 1 00 120 


Time, min 
(c) z axis. 

Figure 3.9.— Typical attitude-rate plot. 
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Figure 3.1 0.— Typical attitude-error plot. 
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Gimballed thruster dynamics 
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4.0 Slosh Model 


One of the major design issues of the COLD-SAT spacecraft was the effect of fluid slosh on spacecraft 
dynamic stability and attitude control system performance because of the large mass of liquid hydrogen 
onboard. Therefore, one of the goals of the COLD-SAT dynamic model was to develop a slosh model 
that would approximate slosh with sufficient accuracy to permit the torques produced on the spacecraft 
by slosh to be evaluated for the conditions that the COLD-SAT spacecraft was expected to encounter. 

An equivalent mechanical model is used to represent the fluid in the experiment tanks. This model 
consists of two independent linear spring-mass-dashpot (SMD) systems per tank: one coupled to the 
spacecraft y axis rotational dynamics and the other to the z axis rotational dynamics. The equivalent 
mechanical model used for slosh is shown in figure 4.1. 


4.1 Computation of Equivalent Mechanical Model Parameters From 
Physical Parameters 

Two orthogonal mechanical SMD systems are used to represent fluid slosh; each is attached to the 
tank wall so that the equilibrium, or null, position of the sloshing mass is on the tank centerline. If 
the proper values of spring constant, mass, damping, and distance below the free surface are selected, 
the mechanical model’s dynamic behavior will under most conditions approximate that of the actual 
sloshing fluid in terms of such externally observable characteristics as reaction force produced on the 
tank wall, slosh frequency, and amplitude. This section is concerned with determining these mechani- 
cal model quantities in terms of actual physical properties, such as tank dimensions, fluid properties, 
fill level, and local acceleration, so that the dynamic behavior of die mechanical model will best 
approximate the actual sloshing fluid (ref. 6). 

Cylindrical geometry is assumed for all tanks. Other assumptions are simplifying assumptions for 
the slosh model: 

(1) Only the fundamental slosh frequency is considered. 

(2) Sloshes in pitch and yaw are assumed independent. 

(3) Bare, smooth tanks are used. 

(4) Slosh can be modeled in any of the three tanks or any combination thereof. 

All parameters are computed at the start of the simulation only and thus are assumed not to change 
significantly during the course of the simulation. 

The equivalent spring constant, K sl (in pounds force per foot), for each tank is computed as follows: 


K sl ~ m T 


1.19/i 


where 


m-i 


if- 


tanh 3.68 _ j + AK sl (a x Ji,d) 


total mass of liquid hydrogen in tank, (l/4>jrp H d 2 h, lbm 


Ph, 


density of liquid hydrogen, lbm/ft 3 
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/ 

d 

A K sl {a x Ji4) 


*tank 


x component of spacecraft linear acceleration produced by axial thrust (F T ) 

x /nl SIC ’ 

ft/sec 2 (For the purposes of the slosh model only, other accelerations in the * 
direction, as well as y-directed and z-directed accelerations, are neglected.) 

nonsloshing height of fluid in tank, x^f, ft 

fill level of tank, percent 
diameter of tank, ft 

correction factor for low acceleration levels, slug/sec 2 or Ibf/ft (For COLD-SAT, 
this term is essentially zero for the three axial-thrust-induced, nonzero acceleration 
levels and is only significant for simulation of slosh under zero-acceleration condi- 
tions. This term cannot be computed analytically and is computed in the model by 
interpolating actual experimental data given in reference 7.) 

height (i.e., x-axis dimension) of tank, ft 


The mass of the equivalent SMD system, m sl (in pounds mass), is computed as follows: 


m sl = (m r)(^j tanh 3 68 ^ 


The distance of the equivalent SMD system below the free surface, k sl (in feet), which causes the 
SMD system dynamic behavior to best approximate that of the fluid, is given by 


K sl 


- t»nh 3.68 A 

3.68 d 


The actual distance of the sloshing mass from the spacecraft CG, D cc (in feet), is given by 


^CG — *tank base 


+ h - L 


where jc^ base is the x coordinate of the tank bottom (in feet). 

For local acceleration levels above zero gravity, the damping ratio of the SMD system, 
(dimensionless), is estimated as follows: 

r rn 4 -^ 8 v°- 5 fA) 075 a °' 25 for bate tank, a > 0 

¥ (2j 1 

where 

2 

v kinematic viscosity of liquid hydrogen, 0.00197 cm /sec 

2 

a x linear acceleration in x direction, cm/sec 
d diameter of tank, cm 
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Note that the value of damping given by this expression becomes infinite as local acceleration 
approaches zero. For values of a x less than 1 |ig, it is assumed that t = 10 -6 . 

The damping coefficient, B si (in pounds force-seconds per foot) for the tank, is used to express 
damping in the slosh dynamical equations given in section 4.3. It is related to the damping ratio as 
follows: 


^sl ^sl ^riySl m sl s i ^ si m si 


where is the slosh natural frequency (in radians per second), 



4.2 Effect of Spacecraft Rotational Motion and Slosh-Model Tank Motion 

The slosh model consists of two independent SMD systems; each system is separately coupled to 
the spacecraft y-axis and z-axis rotational dynamics. Each SMD system is forced by the inertial 
angular displacement (with respect to some arbitrary reference) and the angular rotation about the 
corresponding spacecraft body axis. In turn, the reaction force produced by each SMD system 
produces a torque about the corresponding spacecraft axis. Because the slosh model is a translational 
mechanical system and not a torsional one, it is necessary to develop the relationships between the 
rotational motion of the spacecraft and the corresponding motion of the tank walls, which forces the 
slosh model, as well as the relationship between the reaction forces produced on the tank walls and the 
resulting torque produced on the spacecraft about its CG. The coupling between the spacecraft angular 
motion and the slosh model is given in this section. The coupling between the slosh model and the 
resultant torque produced on the spacecraft is given in section 4.4. 


42.1 Displacement of Tank Wall Due to Spacecraft Rotation 

As the spacecraft rotates about the y or z axis, the null point of the sloshing mass moves along an 
arc of radius D cc . If the angular displacement of the spacecraft is small, this arc can be approximated 
by a straight line tangent to this arc at the null point. Thus, the tank wall is considered to move in a 
straight line as the spacecraft rotates, as given by the following relations: 

v-Axis slosh z-Axis slosh 

>W = D CG®y z tw ~ D CG®z 

where 

Dqq distance of slosh-mass null point from spacecraft CG, ft 

ytw’ z tw y an< * z components, respectively, of linear displacement of null point (and tank walls) 
with respect to inertial space due to spacecraft rotation, ft 

0 y , 0 Z y and z components, respectively, of spacecraft rotation about its CG with respect to 

some arbitrary reference, rad 
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The angles 0 and 0 Z are extracted from the spacecraft attitude quaternion in a manner similar to the 
extraction of Attitude error discussed in section 3.5. The attitude quaternion at the beginning of the 
simulation is used as the attitude reference. The Euler rotation angle extracted from the quaternion 
will always be in the range 0° to 180°, which leads to the restriction that slosh cannot be simulated for 
non-inertially-fixed attitudes: A spacecraft rotation of more than 180°, in any direction, from the 
initial attitude would cause a discontinuity in the computed values of z^, or both, which would 
cause the slosh model to produce incorrect results. The relationship between the Euler rotation axis, 
the Euler rotation angle, and the quaternions is given in appendix B. 

The computation of the angles 0 y and 0 Z from the spacecraft attitude quaternion is given in equation 
form as follows: The slosh rotation quaternion is computed as 

?slosh ~ Q * 


where 

q actual spacecraft attitude quaternion 

< 7 q spacecraft attitude quaternion at start of simulation (reference attitude) 

The slosh rotation angles are given by 

v-Axis slosh z-Axis slosh 

= c 4dos h = ^slosh 

where 

e = <t> = cos' 1 

simp 


4.22 Velocity of Tank Wall Due to Spacecraft Angular Rates 

The slosh SMD system is also forced by the inertial velocity of the tank wall, which is related to 
the spacecraft angular rates about the y and z axes. As was the case for the displacement computation, 
this velocity is actually the tangential velocity at the null point of the sloshing mass. It can be 
computed as follows: 


v-Axis slosh z-Axis slosh 

ytw = *rvv = D CG “z 


where 

ytw’^tw y an ^ ^ components, respectively, of tangential velocity of tank wall, ft/sec 

(0 ,(*) attitude rates about spacecraft y and z axes, respectively, rad/sec 


71 


4.3 Slosh-Model Implementation 


An SMD slosh representation was chosen because of its simplicity and because, unlike a pendulum 
analogy, it is valid at all acceleration levels down to zero. Because of the COLD-SAT spacecraft 
geometry and the orientation of the axial thrust vector, it can validly be assumed that slosh about the 
spacecraft x axis is negligible; thus, the fundamental slosh frequency can be modeled by two indepen- 
dent SMD systems that are coupled to the spacecraft rotational dynamics about the y and z axes, as 
was previously discussed. It is assumed that the slosh harmonics beyond the fundamental are 
insignificant and that there is no coupling between the two SMD systems themselves (although 
coupling exists through the rotation model). Both are valid assumptions if the system is not near 
resonance. 

The mechanical equivalent system, shown inside one tank relative to the fluid free surface, is illustrated 
in figure 4.1. All tanks are assumed to be cylindrical. The parameters m sl , K sl , B sl , and t sl are 

computed from the fluid density, fill level, tank diameter, and linear acceleration (produced by the 
axial thrusters), as was discussed in section 4. 1 . Slosh is assumed to significantly affect the spacecraft 
rotational dynamics only; the effect on the spacecraft trajectory is assumed to be negligible and is not 
fed back to the translation model although lateral motion of the spacecraft is accounted for in the slosh 
model itself. Refer to figure 1 .4 for a block diagram representation of die slosh model. 

The relative displacement of the sloshing mass for tank i with respect to its null point is computed 
as follows: 


v-Axis slosh 


z-Axis slosh 


y rei y tw ~ y si + ys/C z rel ~ z tw z si + Z S/C 


where 

y and z components, respectively, of displacement of sloshing mass with respect to 
inertial space for tank, ft (This is computed by integrating the accelerations resulting 
from forces applied to the sloshing masses.) 

ysiC' z SIC y ant * z components, respectively, of translational displacement of spacecraft due to the 
effects of slosh, ft (This displacement is not fed back to the translation model.) 

The relative velocity of the slosh mass with respect to the tank wall is computed as follows: 

v-Axis slosh z-Axis slosh 

y re i = ytw~ y S i + hie z re i - + *s/c 

where 

y s j>£ S j y and z components, respectively, of velocity of sloshing mass with respect to inertial 

space, ft/sec 

ys/C’Z-sJC y an ^ z components, respectively, of translational velocity of spacecraft due to effects 
of slosh, ft/sec (This velocity is not fed back to the translation model.) 
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The force vector acting on the slosh mass, F st (in pounds force), is computed as follows. The 
negative of this force vector produces a disturbance torque on the spacecraft, as discussed in 
section 4.4, and a translational motion of the spacecraft, which is discussed here. 


?sl = 


K sl}’rel + ^ sl^rel 
^sl z rel + ^sl^rel 


Integrating the acceleration produced by this force yields the velocity and position of this mass with 
respect to inertial space: 


v-Axi8 slosh 

z-Axis slosh 



y s i = & ysftW + *w> 

Z S 1 = tsl&W + Z sl (0) 


The translational acceleration acting on the spacecraft due to slosh in tank i is computed as follows: 


0's/C )« ~ 


v-Axis slosh 

M, 


m S/C ~ ^ 

j 


Os/c )« ~ 


z-Axis slosh 

N, 


m S/C ^ ( m sl)j 

j 


The translational velocities and displacements of the spacecraft due to slosh are computed by 
integrating these accelerations as shown here: 


v-Axis slosh 


z-Axis slosh 


j's/c 


-JC' 


E $s/c¥ T ) 


dj 


L S/C 


-K 


^ C%/c¥ T ) 


idr 


and 


ys/c = £ hlC^ dT z s/c - h* Zs/c( T ) dT 


73 



4.4 Disturbance Torque Produced on Spacecraft by Slosh 


The reaction force on the tank wall F sl was computed in section 4.3. This is the sum of the forces 
produced by the spring and the dashpot; it translates into a torque on the spacecraft. This is the slosh 
torque, which was discussed in section 3.2.2. 


^slosh 52 (^CG )« 

i 


A typical plot of slosh torque (resulting from slosh in all three tanks) for the COLD-SAT spacecraft 
for approximately one orbit, showing x, y, and z components, is given in figure 4.2. 
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(a) x axis. 




(c) z axis. 

Figure 4.2.—Typlcai slosh torque (slosh In all three tanks). 


75 


5.0 Microgravity Accelerations in Experiment Environment 

While on orbit the COLD-SAT spacecraft was to perform a series of experiments to investigate the 
behavior of cryogenic fluids in a near-zero-gravity environment. For some of these experiments it was 
desired to control both the magnitude and direction of the accelerations acting on the spacecraft, 
particularly at the fluid surface, as precisely as possible. In addition, it was desired that the 
acceleration level remain constant throughout the course of the experiment. Three levels of 
acceleration were desired. The acceleration was to be directed along the positive spacecraft x axis, 
with zero acceleration along the y and z axes. Axial thrusters were provided on the aft end of the 
spacecraft for this purpose. The thrusters were sized to provide the following levels of thrust: 0.04, 
0.16, and 0.52 lbf. The acceleration levels generated by these thrust levels will depend upon the 
spacecraft mass at the time that the thrusters are operating. With a spacecraft nominal mass of 
6000 Ibm, the three resulting acceleration levels are 6.7, 26.7, and 86.7 pg directed along the +x axis. 

These desired acceleration conditions are not obtainable in practice because a number of factors act 
to alter the desired acceleration. These undesired accelerations are termed parasitic accelerations and 
are caused by the following conditions: 

(1) Axial thruster misalignment: If the thrust vector generated by the fixed axial thrusters is not 
parallel with the spacecraft x axis, undesired components of acceleration in the y and z directions will 
be produced. 

(2) Gimballed thruster The gimballing action of the gimballed thruster causes the thrust vector to 
change in direction, thereby creating components of thrust in the spacecraft y and z directions and 
causing the x component of thrust to vary in magnitude. 

(3) Spacecraft mass changes: During the course of the mission the spacecraft mass decreases owing 
to consumption of hydrazine and overboard loss of experiment fluids. Thus, for axial thrusters with 
constant thrust level the acceleration level gradually increases throughout the mission. 

(4) Aerodynamic drag: Aerodynamic drag modifies the spacecraft acceleration. The direction of 
this disturbance acceleration, in body coordinates, depends upon the spacecraft attitude selected and the 
spacecraft orbit. 

(5) Spacecraft angular accelerations: Angular accelerations acting on the spacecraft are produced by 
disturbance torques (including slosh) and the action of the attitude control system. Linear 
accelerations result that are tangentially directed with respect to the CG. 

(6) Spacecraft angular rotation: The body rates of the spacecraft produce linear accelerations that 
are radially directed with respect to the CG. 

(7) Gravity gradient: A gravity-gradient acceleration field is produced by the nonuniformity of the 
gravitational field throughout the spacecraft body. 

(8) Uncoupled RCS thrusters: The firing of an uncoupled control thruster, whose function is to 
produce a control torque, will cause linear accelerations of the spacecraft body. For COLD-SAT these 
accelerations are always directed in the +x direction. 
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(9) Angular misalignment of RCS thrusters: If the thrust vectors produced by the RCS are not 
parallel to the principal spacecraft axes because of thruster misalignment, a parasitic linear acceleration 
will be produced in the direction of misalignment when an RCS thruster is fired. 

To gain an understanding of these microgravity disturbances and to evaluate the effect of various 
attitude control systems on the microgravity environment, the rotation model was expanded to include 
the computation of these accelerations at various user-selectable points throughout the spacecraft. A 
qualitative assessment of these accelerations was obtained by plotting the computed acceleration levels, 
typically for one complete orbit. A typical plot of the microgravity environment for a point at the 
front of the supply tank, showing x, y, and z components, is given in figure 5.1 for approximately one 
orbit. The RCS thruster firings are clearly evident, superimposed on the background gravity-gradient 
acceleration. Comparison with figures 3.9, 3.10, and 3.12 reflects a correlation between the 
microgravity environment and the location in the attitude-error/attitude-rate phase plane. 

To obtain a quantitative measure of the degree to which the total accelerations at the selected points 
deviate from the desired acceleration levels, two parameters were defined as follows: 


(1) The maximum instantaneous deviation (MID) 

(2) The root-mean-square deviation (RMSD) 

Precise definitions of these parameters are provided in section 5.5. It was found that these measures 
of miciogravity deviation from the ideal acceleration requirements correlated well with the visual 
impressions created by the graphical representation. For the plot shown in figure 5.1 the maximum 
instantaneous deviation is 19.35 pg along the x axis, 6.916 pg along the y axis, and 8.151 pg along 
the z axis. The RMSD is 1.181 pg along the x axis, 0.7282 pg along the y axis, and 0.2647 pg along 

the z axis. 

In the model the linear accelerations (conditions (1) to (4)) are computed in the translation model, 
and the remaining effects are computed in the rotation model. The acceleration level may be 
computed at as many as 10 points. At each of the 10 points the components of the total acceleration 
along the three body axes is computed. 

The following subsections describe in detail the computation of the various components of the 
microgravity disturbances. Also described are the computation of the MID and RMSD parameters. 


5.1 Microgravity Acceleration Resulting From Spacecraft Linear 
Acceleration 

This component of microgravity acceleration results from translational spacecraft acceleration, 
which is produced by the fixed and gimballed axial thrusters, from aerodynamic drag, and from linear 
forces produced by uncoupled and/or misaligned RCS thrusters: Of these, only the x component of 
axial thrust produces a desired acceleration; the y and z components of axial thrust, as well as all of 
the other causes indicated, produce parasitic accelerations. Computing this acceleration is simply a 
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matter of converting the total external 1 accelerations acting on the spacecraft, which were computed 
by the translation model, into microgravity units. This computation, in equation form, is as follows: 


c = S 


10 ° 


32.175 


f (^xt)s/C> Mg 


where (o^sjc was g' ven 0 n fe® 1 P er second squared) in section 2.5. This acceleration is independent 
of location within the spacecraft. 


5.2 Microgravity Acceleration Resulting From Spacecraft Rotational 
Dynamics and Gravity-Gradient Effects 

The remaining three components of microgravity acceleration are unique in the sense that the 
magnitudes and directions of these accelerations at a given point in the spacecraft are a function of the 
coordinates of that point with respect to the CG. Therefore, it is necessary to define the location at 
which the microgravity accelerations are to be computed. This is done by means of a vector from the 

CG to the point of microgravity computation, denoted r* p , which is used in all of the acceleration 
computations that follow. The microgravity accelerations discussed in this section are entirely 
parasitic. 


5.2.1 Angular-Rate-Induced Centrifugal Acceleration 

Owing to centrifugal effects, spacecraft rotation about the CG produces a radially directed 
acceleration that is a function of angular rates and relative distance from the CG. Because 
COLD-SAT’s attitude control scheme results in extremely low angular rates, this acceleration is 
typically insignificant. However, it is included in the model for generality. 

The centrifugal microgravity acceleration is computed as follows: 


cen 


io b 

132.175 


f C ? x (u? x rp, pg 


where 

u? spacecraft angular rate vector, rad/sec 

vector from CG to point of computation, ft 


^Because the spacecraft is in free fall, the gravitational acceleration does not affect the microgravity environment and thus is not 
included in the computation. 
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5.2 2 Angular-Acceleration-Induced Tangential Acceleration 

Spacecraft angular acceleration about its CG produces a tangentially directed linear acceleration that 
is a function of the distance from the CG and the spacecraft angular acceleration vector. This 
acceleration is produced by all torques that act on the spacecraft, including both disturbance and 
control torques. For COLD— SAT, because of the large magnitude of the RCS control torques with 
respect to other torques acting on the spacecraft, this component of microgravity acceleration is 
dominant when RCS thrusters are firing. This acceleration is predominantly in the form of large- 
magnitude, short-duration pulses that correspond to RCS activity. This effect can be clearly seen in 
figure 5.1. 

The tangential microgravity acceleration is computed as follows: 


J tan 


10 b 

132.175 


f w* r}„ pg 


where to* the spacecraft angular acceleration vector (in radians per second squared). 


5.23 Microgravity Acceleration Resulting From Gravity-Gradient Effects 

Gravity-gradient microgravity acceleration results from the fact that (like the gravity-gradient 
disturbance torque acting on the spacecraft rigid body) only the CG of die spacecraft is truly in free 
fall and nowhere else inside the spacecraft does the orbital centrifugal acceleration (away from the 
Earth) exactly balance the gravitational acceleration (toward the Earth). This results in a “gravity- 
gradient field” that is zero at the CG and increases as one moves away from the CG. 2 The variation 
of the field is different along the direction of the spacecraft position vector R than it is in the 
remaining two orthogonal coordinate directions. 


Computation of the gravity-gradient acceleration requires that the vector r* p be transformed into the 
local-horizon (LH) coordinate system. This system has its origin located at the spacecraft CG, with 
the z axis in the direction of the zenith and the y axis perpendicular to the spacecraft orbit plane. The 
x axis is chosen to form a right-handed system and will coincide with the velocity vector for a circular 
orbit only. 

The vector r* p in spacecraft body coordinates is transformed into TOD coordinates (fppoD ^y 
means of vector transformation and by using the spacecraft actual attitude quaternion as follows: 


0 


0 

1 

V 

8 

= 9* 

fp . 


where q is the attitude quaternion. 


2 Note that the microgravity acceleration that is produced by the gravity-gradient field is in addition to the microgravity acceleration that 
is caused by gravity-gradient-disturbance-torque-induced angular acceleration on the spacecraft, as previously discussed. 
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The transformation between TOD coordinates and LH coordinates is now derived. The LH system 
can easily be defined in terms of previously defined base vectors (see section 2.6.3) as follows: The 
LH z axis is the unit position vector v* R . The LH y axis is a unit vector normal to the orbit plane, 
v* N = v* R * iT v , where is the unit velocity vector. The LH x axis is the crossproduct * v* R . 
Thus, the transformation matrix from TOD to LH coordinates is 

[TOD — * LH] = * v* R | v* N | ujjf] 


The vector f* p can now be expressed in LH coordinates (rp LH as 

<V»LH * [TOD - LHK/^W 

The microgravity acceleration vector due to gravity-gradient effects, denoted i^ G , is computed in LH 
coordinates as follows: 


( *gg)lh " 


10 e 


32.175 


Earth 


R- 


px 


py 


-2 r 


pz. 


> Mg 


LH 


where 

R magnitude of spacecraft position vector, ft 

MEaith Earth gravitational constant, 1 .4076469* 10 16 frVsec 2 


The acceleration vector is now transformed to TOD coordinates (sgg>tod an< * finally to spacecraft 
body coordinates (i^) as follows: 

WfG^TOD = [TOD — > LHJ^Sqq)^ 


0 

= (q c )* 

0 



*GG 


( s GGhoD 


5.3 Total Microgravity Acceleration 

The total microgravity acceleration is the vector sum of the four components: 

r = K- + *GG + £en + C 
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5.4 Plotting Total Acceleration 


It is required to graphically plot the total microgravity acceleration, which is the sum of the 
translational, gravity-gradient, centrifugal, and tangential accelerations as previously discussed. When 
there is no RCS activity, the dominant torque acting on the spacecraft is the gravity-gradient torque, 
and perhaps the reaction-wheel- and gimballed-thruster-produced torques if these systems are included 
in the simulation. In any event the torques produced by these three entities are smooth curves that can 
easily be plotted. However, if one or more RCS thrusters are fired, the torques produced by these 
thrusters are by far the largest torques acting on the spacecraft for the duration of the thruster on-time. 
But because this on-time is very short, typically with a value of t min on (see section 3.2.3.2) or slightly 
greater, the resulting tangential acceleration profile consists of very narrow pulses of large magnitude 
whose width matches the thruster on-time. These pulses are too short in duration to plot in most cases 
(typical spacing between plot points is on the order of 5 to 10 sec, whereas on for COLD-SAT is 
200 msec). In theory, the number of data points plotted could be increased so that the spacing 
between points is on the order of 200 msec, but this would result in excessively large data files. A 
better solution to this problem is to exaggerate the width of the RCS pulses (for plotting purposes 
only) from their actual width of about 200 msec to several seconds so that the resulting plots will 
faithfully illustrate the microgravity environment. 

Because these short pulses occur only when an RCS thruster is fired, the algorithm implemented 
simply checks the states of all RCS thrusters; when one or more RCS thruster is in a full-on state (i.e., 
it has passed the startup transient), the algorithm holds the value of microgravity acceleration 
computed at this point (which is die maximum value) for time r hoId , which may be several seconds, 
and allows the plotting program to catch at least one point at this maximum value. After /j,old seconds 
have passed, the value of the stretched acceleration variable is again set to that computed by the 
simulation. This results in pulses that are clearly visible on the plotted output. The value of ^old ' s 
selected so that r hold > t mc , where is the time between plotted points. 


5.5 Computation of Microgravity Deviation Parameters 

It is desirable to develop a numerical measure of how much the actual microgravity environment 
deviates from the ideal microgravity acceleration for each point. Two vectors are defined for each 
point where the microgravity environment is to be computed: 

(1) Maximum instantaneous deviation (MID): This parameter represents the largest magnitude of 
deviation of the actual acceleration from the desired acceleration for one orbit. 


(2) Root-mean-square deviation (RMSD): This parameter represents an average deviation of the 
actual acceleration from the desired value over one orbit. 

The MID and RMS deviations are expressed as follows: 
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where 

s nom 
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RMSD 
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[5^(1) •tnonJ 
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' <h 
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V 

S*( T) 
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Mg, over one orbit 


nominal acceleration in +x direction (computed at beginning of simulation only), 
nom /( m S/c) nom 

nominal thrust level (0.04, 0.16, or 0.52 lbf) in +jc direction 
nominal spacecraft mass, slug 


12 

8 

4 

0 



(a) x axis. 



Time, min 


(c) z axis. 


Figure 5.1 .—Typical microgravity plot. (Point of computation 
Is located on x axis at front of supply tank.) 
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Appendix A 
Notation and Constants 


TABLE A. 1— NOTATION CONVENTIONS 


Symbol 

Typical usage 


unit vector; subscript indicates which one (in this case, v subscript indicates 
unit velocity vector, R subscript indicates unit position vector, etc. 

*<, 

unit vectors along coordinate axis as indicated by subscript 

v, R ,a t etc. 

vector (in general, not unit vector) 

v r V v z 

x t y, and z components of vector V 

7 r. V 7 e 

R , and 0 components of vector V (spherical coordinates) 

^S/C 

vector referenced to coordinate system other than one in which vector was 
originally defined (e.g., is spacecraft velocity vector expressed in 

spacecraft coordinates; original definition of this vector was in TOD 
coordinates) 

$ 

quaternion 

,<0) 

scalar part of quaternion # 

«». <?<» 4® 

vector part (*, y, and z components, respectively) of quaternion $ 

0 

three-axis attitude error (in strict mathematical sense, not a vector) 

I 

rectangular matrix 

'ij 

single element of matrix I (/th row, jth column) 

I T 

transpose of matrix I 

r* 

inverse of square matrix I 

x. If. ^ 

time derivative (rate of change) of scalar, vector, and quaternion, 
respectively 

[A - B1 

transformation matrix from coordinate system A to coordinate system B 
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TABLE A.2.— FUNDAMENTAL PHYSICAL CONSTANTS 


Symbol 

Name 

Value 

B 

Earth gravitational constant (product 
of universal gravitational constant 
and Earth mass) 

1.4076469* 10 16 f^/sec 2 
3.986005* 10 14 m 3 /sec J 


Earth gravitational second zonal 
harmonic coefficient 

1.0827* 10‘ 3 

h 

Earth gravitational third zonal 
harmonic coefficient 

-2.56* HT 6 

j* 

Earth gravitational fourth zonal 
harmonic coefficient 

-1.58* 10" 4 

*Earth 

Earth equatorial radius 

20 925 656 ft 
3443.920 n mi 
6378.15 km 

5 Sun 

Solar radius (visible surface) 

2.283436*10? ft 

^Emih-Svm 

Earth -Sun distance (semimajor axis), 
AU 

4.908066* 10* 1 ft 

M E**h 

Earth magnetic dipole strength 

8.1 MO 22 A-m J 

MO 

Magnetic permeability 

4ir* 10 -7 H/m 


84 










Appendix B 
Quaternion Algebra 

Because quaternions are used throughout this model to represent spacecraft attitude, it is essential to 
provide some background information on the mathematics of quaternions. 


B.l Introduction 

A quaternion is an ordered quadruple that appears to be a complex number with an imaginary part 
consisting of a three-element vector. Quaternions can be expressed in many different forms, but two 
forms are of concern here. The unit-vector form is expressed as follows. 




,<‘>r 




,<»r 


where CyT and j T are unit vectors similar to those encountered in vector analysis. Alternatively, a 
quaternion can be expressed in column-vector form as 


9 


9 

9 

9 

9 


<o> 

0) 

< 2 > 

< 3 > 


where q 0 is called the scalar part and q v q 2 , and q 3 make up the vector part. The column-vector form 
is used throughout this document when referring to quaternions in general. However, the unit-vector 
form is convenient when dealing with quaternion multiplication, as will be shown shortly. 

A special type of quaternion is the unit quaternion. The unit quaternion simply satisfies the 
property that 


? <2 f * (* 0> F - ■ 


Because many of the properties defined herein are only valid for unit quaternions, all quaternions used 
in this model, unless otherwise indicated, are assumed to be unit quaternions. 


B.2 Advantages of Quaternion Attitude Representation 

A unit quaternion actually represents a single three-dimensional rotation of a given angle about a 
given vector. Thus, a unit quaternion can be used to express a transformation between two different 
coordinate systems and is used to represent the attitude of the spacecraft in the COLD— SAT dynamic 
model. 
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Alternative representations of spacecraft attitude are the use of Euler angles and the attitude matrix 
(sometimes called a direction-cosine matrix). Euler angles represent three sequential two-dimensional 
rotations to describe a rotation in three dimensions. The equations of motion for the Euler angle 
representation contain nonlinear functions. Singularities occur for certain rotation angles, making the 
Euler angle approach valid only for small-angle rotations. The attitude matrix, on the other hand, has 
linear equations of motion but requires integration of nine values per time step, as opposed to three for 
the Euler angle representation. Another problem with the attitude matrix representation is the tendency 
of the matrix to drift away from orthonormality because of numerical inaccuracies in the integration. 
Restoration of orthonormality is a nontrivial procedure that requires many operations. The quaternion 
representation suffers from none of the aforementioned problems. It requires only four integrations per 
time step, is valid for any rotation, and has no singularities anywhere in the equations of motion. The 
only condition that must be maintained is normality of the quaternion, which is easily checked and 
corrected. 


B.3 Quaternion Mathematics 

B3.1 Addition 


Quaternions are added together just like vectors or column matrices. If p and q are quaternions, 
the quaternion sum is found as follows: 


> 


> 


■„< 0> 

+ q w> 

P <1> 


«® 


P<1) 



T 

«® 


p® 

♦ f <*> 

p°\ 


s 0> 


p 0> 

+ q Q> 


B3.2 Multiplication by Scalar 

Quaternions can be multiplied by a scalar like vectors and matrices. If q is a quaternion and a is a 
scalar, then 




aq*® 

»® 


aq ^ 

»® 


aq ® 

jr® 


aq » 
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B33 Quaternion Multiplication 

If p and q are quaternions defined as 


,< 0 ) 

,<*> 

,< 2 > 

,<3> 


= pQb + p^k 


and 


= q<o> + ,<»r+ q v>r* <7 0> ** 


where i ? h and ** are unit vectors, multiplication of the quaternions p and q, denoted p * q , is 
defined as follows by using the unit-vector representation. 

(1) Multiply each element of p by each element of q, preserving the order of multiplication 
(16 multiplications). 

(2) Simplify the resulting terms by using the following rules. Note that this operation is not 
commutative: 

nr= -i j*j* = -i k*k = - i 

nr = k* = r **/ = j 

ic*f = -V j*i* = ~k* i*ic* = -j 


<7 = 


,<o> 

,< 1 > 

I® 

, 0 > 


Note the similarity to complex number multiplication and vector crossproduct. 


B3.4 Conjugation 

If q is a quaternion, the conjugate of q, denoted by q c , is found by negating the vector part of the 
quaternion q: 
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"<7 

-9 

."9 


<l> 

< 2 > 

<3> 


Note the similarity to the conjugation of a complex number. 


B3.5 Inverse 


The inverse of a quaternion q, another quaternion denoted by q 1 , satisfies the following property 




= 1 


If q is a unit quaternion, then q 1 * q c (i.e., the inverse and conjugate of any unit quaternion are 
equal. 


B3.6 Extraction of Rotation Angle and Rotation Vector 

Because a quaternion represents a rotation of a given angle about a given vector, the rotation angle 
6 and rotation vector 0 - [v^ v v z ] that a given unit quaternion q represents can be depicted as 
follows: 


9 = 


r <°> 

r°> 

r® 

,<3> 


0 


cos. 


. e 

v_ sm_ 
* 2 

. 0 
v y sm j 

. 0 

v. sin— 

z 2 


The rotation angle 0 is found from the scalar part of the quaternion. If 0 is restricted to two quad- 
rants, there is a one-to-one correspondence between the quaternion and the rotation in space. If the 
first and second quadrants (0 £ 0 £ n) are arbitrarily chosen, then 0 can be solved for. The rotation 
vector is then easily determined from the vector part of the quaternion. 
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B 3.7 Vector-Coordinate Transformations Using Quaternions 

nr 

jf y*= [v v v is a vector defined in TOD coordinates, w = [w x w^\ is the same 

vector expressed in spacecraft body coordinates, and q = [<7^ q^ q^ q^ * s the spacecraft 
attitude quaternion, the vector-coordinate transformation from TOD to spacecraft body coordinates can 
be expressed as 



This expression in compact form is 



The inverse transformation from spacecraft body to TOD coordinates is 
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Appendix C 
Coordinate Systems 

C.l TOD System 

The TOD (geocentric equatorial true equinox of date) coordinate system (denoted X, Y, and Z) is 
a pseudo-inertial 3 right-handed Cartesian system with its origin situated at the center of the Earth. 

The X-Y plane is coincident with the Earth’s equatorial plane and the Z axis passes through the North 
Pole. The X axis points in the direction of the vernal equinox (first point of Aries) on the date in 
question. This coordinate system is illustrated in figure C.l. Vectors expressed in this system are 
given the subscript TOD. This coordinate system is the basis for the translation model. 


C.2 Spacecraft Body System 

The spacecraft body coordinate system (denoted x, y, and z) is a right-handed Cartesian coordinate 
system whose origin is defined as the instantaneous, actual CG point of the spacecraft. The x axis is 
parallel to the long axis of the spacecraft, the positive direction being in the direction of the small 
receiver tank. The y axis is defined to be parallel to the solar panel connecting struts, and the z axis is 
in the direction opposite to the parabolic antenna. This coordinate system is not fixed with respect to 
the spacecraft rigid body but moves with the CG location. This coordinate system is the basis for the 
rotation model because all attitudinal motion of the spacecraft will occur about the CG. This system is 
illustrated in figure C.2. Vectors expressed in this system are given the subscript SAD. 


C.3 Alternative Spacecraft Body System 

The alternative spacecraft body coordinate system (denoted jc', y', and z!) is a right-handed 
coordinate system similar to the spacecraft body system discussed previously. The x!, /, and z' axes 
point in the same directions, respectively, as the x, y, and z axes, but the origin of this coordinate 
system is located at the intersection of the plane containing the back thrusters and the line perpen- 
dicular to this plane that passes through the nominal spacecraft CG. Note that this system, unlike the 
x, y, z system, is fixed with respect to the spacecraft rigid body. It is thus convenient to use it to 
specify the locations of thrusters, tanks, and other components that are fixed with respect to the 
spacecraft rigid body. This system is illustrated in figure C.2. 


C.4 Spacecraft Orbit Plane System 

The spacecraft orbit plane (SCO) system is used in the initialization routine that converts a set of 
user-specified orbital elements into an equivalent position and velocity vector. This system is briefly 
mentioned in section 2.2 and defined as an Earth-centered, right-handed Cartesian coordinate system 
that is fixed with respect to the perigee of the spacecraft orbit. The x-y plane of this system is 
coincident with the spacecraft orbit plane, with the x axis passing through perigee of the orbit (if the 


Pseudo-inertial refers to the fact that this coordinate system does move very slowly with respect to inertial space because of the 
precession and nutation of the Earth’s axis. However, over the proposed life of the COLD-SAT spacecraft, this movement is so slight as to 
be negligible. Therefore, for these purposes the TOD system can be considered to be an inertial frame. 
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orbit is circular, the orientation of the x axis is arbitrary). The z axis is in the direction of the space- 
craft angular momentum vector. Vectors expressed in this system are given the subscript SCO. 


C.5 Earth Magnetic Field System 

The Earth magnetic field coordinate system (denoted X,YZ) is a right-handed Cartesian coordinate 
system that is fixed with respect to the rotating Earth. The z axis of this system is defined as passing 
through the Earth’s north magnetic pole, and the x-y plane contains the magnetic equator, with the 
x axis passing through the meridian of longitude in which the north magnetic pole lies. Vectors 
expressed in this system are denoted with the subscript M. This system is used in computing the 
magnetic torque. 

The transformation from the true-of-date (TOD) system to the magnetic (M) system consists of two 
successive Euler rotations. Starting with the TOD system, the first rotation is about the Z axis (i.e., 
the Earth’s axis) by an angle 5, where 6 is the angle between the vernal equinox (X axis) and a vector 
from the Earth’s center through the point on the Earth’s surface where the magnetic equator intersects 
the meridian of longitude containing the north magnetic pole. This intermediate system is denoted Af. 
The angle B is defined as follows: 


B = a - GHA 


where a is the longitude west of Greenwich of the north magnetic pole and GHA is the Greenwich 
hour angle, the angle measured westward from the Greenwich meridian to the vernal equinox. The 
Euler rotation matrix for this first rotation is 


[TOD - A/'] 


cos B -sin 8 0 
sin B cos B 0 
0 0 1 


For the second Euler rotation the A f system is rotated about its y axis by angle P to form the M 
system, where P is defined as 90° minus the latitude of the north magnetic pole. This second 
transformation can be expressed as 


[Af - M] 


COS P 

0 

-sin P 

0 

1 

0 

sin p 

0 

COS P 


Therefore, the complete transformation from TOD coordinates to magnetic coordinates is given by the 
following matrix: 
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cos P cos 5 -cos P sin 6 -sin P 

[TOD -* M] = sin 8 cos 6 0 

sinP cos 5 -sin P sin 5 cos P 

The values of a and P vary somewhat over time as the Earth’s magnetic poles move but are generally 

known for a given date. The GHA is a simple function of the date and time: 

GHA = GHAq + R^t 

where 

GHAq Greenwich hour angle at start of simulation (found in Astronomical Almanac) 

^Eaith Earth rotation rate, deg/sec 
t simulation time, sec 

C.6 Ecliptic Coordinate System 

The ecliptic coordinate system is a heliocentric, right-handed Cartesian coordinate system whose 
x-y plane is coincident with the plane of the ecliptic, whose x axis points in the direction of the first 
point of Aries, and whose z axis is in the direction of the angular momentum vector of the Earth as it 
revolves about the Sun. Vectors expressed in this system are given the subscript e. This system is 
used in computing the desired spacecraft attitude for certain attitude options. 


z 

t ,-North Pole 



X (vernal equinox) 

Figure C.1 . — TOD coordinate system. 
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— Solar panel 


Figure C.2. —Spacecraft body coordinate system. 






Solar panel 
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Appendix D 

Implementation of Model 

D.l Overview 

The COLD-SAT dynamic model is implemented by using the EASY5 engineering analysis system. 
EASY5 is an analysis and design tool for simulation and control systems that is published by Boeing 
Computer Services. It incorporates the following features: 

(1) The convenient block-diagram-based system modeling approach subdivides the model into 
separate components that can be developed and tested independently. 

(2) The transient solution of systems is described by linear or nonlinear ordinary differential 
equations in time. One of several fixed-step or variable-step numerical integration schemes is used, 
including the Runge-Kutta, modified stiff-gear, Adams-Bashforth predictor/Adams-Moulton corrector, 
and Euler methods. The variable-step Adams-Bashforth/Adams-Moulton method was used for all 
COLD-SAT dynamic model simulations, for it yielded excellent stability and accuracy while 
consuming the least CPU time. 


(3) “Standard component” libraries are available that contain a rich assortment of common control 
system elements, including basic first- and second-order transfer functions, nonlinear effects, basic 
control laws, and tabular function generators. Furthermore, user-specified Fortran and Macro 
components can be created to model additional effects by inclusion of the appropriate Fortran code. 
“Fortran components” were used to model the vast majority of COLD-SAT dynamic model 
subsystems. 

(4) Linear control system analysis and design tools, such as root locus, frequency response, and 
linear-model generation, can be performed on both linear and nonlinear models. (Nonlinear models 
are automatically linearized about a given operating point.) 

(5) High-resolution graphical display and hard copy of analysis may be output. 

(6) Model code can call external subroutines. This approach, which permits modular program 
development, was extensively employed for the COLD-SAT dynamic model. 

All of these capabilities are integrated into a single package that runs on a variety of platforms, 
including the Apollo Domain workstations, the VAX/VMS minicomputer, and the Cray X/MP 
supercomputer. Figure D.l is a pictorial overview of the EASY5 system. For additional information, 
see reference 8. 


D.2 Subroutines 

The COLD-SAT dynamic model utilizes several external subroutines and functions. Some were 
developed specifically for this model, but others are general-purpose routines. 
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D 2.1 COLD-SAT Model-Specific Subroutines 


SCATT (ATTUD) computes desired spacecraft attitude direction-cosine matrix for given attitude 

option, axis alignment, optional rotation, and time (section 2.6). ATTUD is 
called by SCATT. 

SLOSHP computes equivalent mechanical slosh-model parameters for liquid hydrogen, 

by using an SMD slosh model, from tank dimensions, fill level, and local 
acceleration (section 4.1) 

TTLACC computes microgravity acceleration for given point in the spacecraft 

(section 5.0) 

IRCSTGL function that returns nonzero value when any RCS thruster is fired. This 

value is used in the microgravity acceleration plotting routine (section 5.4). 

D2J General Astronomical Routines 

ATMJ (FUNC1, function that returns atmospheric density at given point in spacecraft orbit by 

FUNC, and SIMPK) using Jacchia’s 1970 model (ref. 2). FUNC1, FUNC, and SIMPK are called 

by ATMJ. 

JDATE computes Julian date with respect to year 2000 from given month, day, year, 

and Greenwich mean time 

ORBEL computes classical elliptized orbital elements from position and velocity 

vectors 

SGHA computes Greenwich hour angle for given date and GMT 

SPHGGT computes gravity-gradient torque by using spherical Earth model 

(section 3.2. 1.1) 

STVEC2 computes position and velocity vectors from given orbital elements 

(section 2.2) 

SUNVEC computes position of Sun in TOD coordinates on given Julian date, as well 

as vector normal to ecliptic (section 2.6.2) 

SUNVIS determines fraction of Sun’s disk that is visible for given point in spacecraft 

orbit (section 3.2. 1.4) 


D2J General Vector, Matrix, and Other Mathematical Routines 

AROT computes matrix that performs Euler rotation of angle 6 about given axis 

MPROD1 computes product of 3 x 3 matrix and a three-element vector: v 2 = Avj 
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MPROD2 

computes product of transpose of 3 * 3 matrix and a three-element vector 
v 2 = (AVi 

MTMPY1 

computes product of two matrices: C nxp = A nXOT B TOXp 

MTMPY2 

T 

computes product of two matrices: C nxp = A nxm (H pxm ) 

MTMPY3 

T 

computes product of two matrices: C„ xp = (A,,,,,,,) B mxp 

SP2REC 

converts vector in spherical coordinates to Cartesian coordinates, provided 
that <(> component (longitudinal component) of vector is zero 

VCRS 

computes crossproduct of two vectors in three dimensions: w = v * u 

VDOT 

computes dot product of two vectors: x = v • u 

VTSCAL 

computes product of vector and scalar, w = av 

XUNIT 

computes magnitude of given vector, as well as unit vector in same direction 
as given vector 

HYPOT 

/ 2 2 

computes hypotenuse of right triangle: z = \x + y 

D2.4 Quaternion Routines 

DCMTOQ 

converts direction-cosine (attitude) matrix into equivalent attitude quaternion 
(section 2.6.8) 

Q2RPY 

extracts rotation angles about spacecraft x, y, and z axes from quaternion 
specifying rotation from some reference attitude (sections 3.5 and 4.2.1) 

QCVQ 

converts vector from coordinate system A to coordinate system B by using 
quaternion (appendix B, section B.3.7) 

QVQC 

converts vector from coordinate system B to coordinate system A by using 
quaternion (appendix B, section B.3.7) 

QMULT2 

computes product of two quaternions: p = q x *q^ (appendix B, 
section B.3.3) 

QMULT3 

computes product of two quaternions: p = (appendix B, 

section B.3.3) 
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D.3 Code Listing 

D.3.1 EASY5 Model 

The COLD-SAT dynamic model was developed and debugged on an Apollo workstation by using 
the EASY5/W program. EASY5/W is a graphical preprocessor to EASY5 that allows the user to 
create and connect both standard and user-defined components in a manner similar to a computer- 
assisted design (CAD) system. Once a valid model is created in EASY5/W, a batch-format EASY5 
model-generation source file is created from the graphical model. This file can be compiled and run 
on the Apollo or ported to another machine (with the proper connectivity), such as the Cray, for faster 
execution. The output files can then be sent back to the Apollo or to a personal computer for 
graphical postprocessing (see figure D.l). 

Figure D.2 is the EASY5/W Apollo workstation block diagram of the entire model. 


DJ2 Subroutines 

In order to modularize the code as much as possible, much of the code used to implement the 
model was written in the form of standard Fortran 77 subroutines external to EASY5. These 
subroutines are called by the main model. 



Figure D.l . — EASY5 modeling system. 
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Figure D.2. — Apollo workstation block 
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diagram of COLD-SAT dynamic model. 
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Appendix E 

Acronyms and Abbreviations 

A ampere (a unit of electric current) 

AU astronomical unit 

ACS attitude control system 

CG center of gravity of spacecraft 

CPU central processing unit 

COLD-SAT Cryogenic On-orbit Liquid Depot — Storage, Acquisition, and Transfer 
pg unit of acceleration (32.16 x 10~ 6 ft/sec 2 ) 

GHA Greenwich hour angle 

GMT Greenwich mean time 

H henry (a unit of inductance) 

lbf pound force 

lbm pound mass 

MID maximum instantaneous deviation 

m meter (a unit of length) 

N newton (a unit of force) 

N 2 H 4 hydrazine 

P-I proportional-integral (control law) 

psia pounds per square inch absolute 

RCS reaction control system 

RMSD root-mean-square deviation 

S/C spacecraft (also refers to primary (x,yx) spacecraft coordinate system) 

slug 32.175 lbm 

SCO spacecraft orbit system 
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SMD 

TOD 

TOL 

V 


spring-mass-dashpot (slosh) model 
true-of-date coordinate system 
tolerance 

volt (a unit of electric potential) 
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